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:

Public Member Functions
	DisContField ()
	Default constructor.
	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.
	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.
	DisContField (const DisContField &In, const bool DeclareCoeffPhysArrays=true)
	Constructs a 1D discontinuous field based on an existing field.
	DisContField (const DisContField &In, const SpatialDomains::MeshGraphSharedPtr &graph, const std::string &variable, const bool SetUpJustDG=true, 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.
	~DisContField () override
	Destructor.
GlobalLinSysSharedPtr	GetGlobalBndLinSys (const GlobalLinSysKey &mkey)
	For a given key, returns the associated global linear system.
bool	SameTypeOfBoundaryConditions (const DisContField &In)
	Check to see if expansion has the same BCs as In.
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} \).
void	EvaluateHDGPostProcessing (const Array< OneD, const NekDouble > &coeffs, Array< OneD, NekDouble > &outarray)
	Evaluate HDG post-processing to increase polynomial order of solution.
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.
unsigned	GetTraceElmtId (const unsigned elmtid, const unsigned traceid)
ExpListSharedPtr &	GetLocElmtTrace ()
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.
Public Member Functions inherited from Nektar::MultiRegions::ExpList
	ExpList (const ExpansionType Type=eNoType)
	The default constructor using a type.
	ExpList (const ExpList &in, const bool DeclareCoeffPhysArrays=true)
	The copy constructor.
	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.
	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.
	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.
	ExpList (SpatialDomains::PointGeom *geom)
	Specialised constructors for 0D Expansions Wrapper around LocalRegion::PointExp - used in PrePacing.cpp.
	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.
	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.
	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.
virtual	~ExpList ()
	The default destructor.
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\).
int	GetNcoeffs (const int eid) const
	Returns the total number of local degrees of freedom for element eid.
ExpansionType	GetExpType (void)
	Returns the type of the expansion.
void	SetExpType (ExpansionType Type)
	Returns the type of the expansion.
int	EvalBasisNumModesMax (void) const
	Evaulates the maximum number of modes in the elemental basis order over all elements.
const Array< OneD, int >	EvalBasisNumModesMaxPerExp (void) const
	Returns the vector of the number of modes in the elemental basis order over all elements.
int	GetTotPoints (void) const
	Returns the total number of quadrature points m_npoints \(=Q_{\mathrm{tot}}\).
int	GetTotPoints (const int eid) const
	Returns the total number of quadrature points for eid's element \(=Q_{\mathrm{tot}}\).
int	GetNpoints (void) const
	Returns the total number of quadrature points m_npoints \(=Q_{\mathrm{tot}}\).
int	Get1DScaledTotPoints (const NekDouble scale) const
	Returns the total number of qudature points scaled by the factor scale on each 1D direction.
void	SetWaveSpace (const bool wavespace)
	Sets the wave space to the one of the possible configuration true or false.
void	SetModifiedBasis (const bool modbasis)
	Set Modified Basis for the stability analysis.
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.
void	SetPhys (int i, NekDouble val)
	Set the i th value of m_phys to value val.
void	SetPhys (const Array< OneD, const NekDouble > &inarray)
	Fills the array m_phys.
void	SetPhysArray (Array< OneD, NekDouble > &inarray)
	Sets the array m_phys.
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.
bool	GetPhysState (void) const
	This function indicates whether the array of physical values \(\boldsymbol{u}_l\) (implemented as m_phys) is filled or not.
void	MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
	multiply the metric jacobi and quadrature weights
void	DivideByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
	Divided by the metric jacobi and quadrature weights.
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})\).
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.
void	IProductWRTDirectionalDerivBase (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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})\).
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.
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.
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.
GlobalLinSysKey	HelmSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff=StdRegions::NullVarCoeffMap, const StdRegions::VarFactorsMap &varfactors=StdRegions::NullVarFactorsMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray, const bool PhysSpaceForcing=true)
	Solve helmholtz problem.
GlobalLinSysKey	LinearAdvectionDiffusionReactionSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff=StdRegions::NullVarCoeffMap, const StdRegions::VarFactorsMap &varfactors=StdRegions::NullVarFactorsMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray, const bool PhysSpaceForcing=true)
	Solve Advection Diffusion Reaction.
GlobalLinSysKey	LinearAdvectionReactionSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff=StdRegions::NullVarCoeffMap, const StdRegions::VarFactorsMap &varfactors=StdRegions::NullVarFactorsMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray, const bool PhysSpaceForcing=true)
	Solve Advection Diffusion Reaction.
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.
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\).
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.
void	Reset ()
	Reset geometry information and reset matrices.
void	ResetMatrices ()
	Reset matrices.
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.
void	SetCoeff (int i, NekDouble val)
	Set the i th coefficiient in m_coeffs to value val.
void	SetCoeffs (int i, NekDouble val)
	Set the i th coefficiient in m_coeffs to value val.
void	SetCoeffsArray (Array< OneD, NekDouble > &inarray)
	Set the m_coeffs array to inarray.
int	GetShapeDimension ()
	This function returns the dimension of the shape of the element eid.
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.
void	ImposeDirichletConditions (Array< OneD, NekDouble > &outarray)
	Impose Dirichlet Boundary Conditions onto Array.
void	ImposeNeumannConditions (Array< OneD, NekDouble > &outarray)
	Add Neumann Boundary Condition forcing to Array.
void	ImposeRobinConditions (Array< OneD, NekDouble > &outarray)
	Add Robin Boundary Condition forcing to Array.
void	FillBndCondFromField (const Array< OneD, NekDouble > coeffs)
	Fill Bnd Condition expansion from the values stored in expansion.
void	FillBndCondFromField (const int nreg, const Array< OneD, NekDouble > coeffs)
	Fill Bnd Condition expansion in nreg from the values stored in expansion.
void	AvgAssemble (bool useComm=true)
	Assemble the average global coefficients \(\boldsymbol{\hat{u}}_g\) from the local coefficients \(\boldsymbol{\hat{u}}_l\) .
void	AvgAssemble (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool useComm=true)
void	LocalToGlobal (bool useComm=true)
	Gathers the global coefficients \(\boldsymbol{\hat{u}}_g\) from the local coefficients \(\boldsymbol{\hat{u}}_l\).
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\).
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.
NekDouble	GetCoeffs (int i)
	Get the i th value (coefficient) of m_coeffs.
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.
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.
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.
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.
NekDouble	Integral ()
	Calculates the \(H^1\) error of the global spectral/hp element approximation.
NekDouble	Integral (const Array< OneD, const NekDouble > &inarray)
NekDouble	VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &inarray)
Array< OneD, const NekDouble >	HomogeneousEnergy (void)
	This function calculates the energy associated with each one of the modesof a 3D homogeneous nD expansion.
void	SetHomo1DSpecVanVisc (Array< OneD, NekDouble > visc)
	This function sets the Spectral Vanishing Viscosity in homogeneous1D expansion.
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.
LibUtilities::TranspositionSharedPtr	GetTransposition (void)
	This function returns the transposition class associated with the homogeneous expansion.
NekDouble	GetHomoLen (void)
	This function returns the Width of homogeneous direction associated with the homogeneous expansion.
void	SetHomoLen (const NekDouble lhom)
	This function sets the Width of homogeneous direction associated with the homogeneous expansion.
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.
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.
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.
int	GetExpSize (void)
	This function returns the number of elements in the expansion.
size_t	GetNumElmts (void)
	This function returns the number of elements in the expansion which may be different for a homogeoenous extended expansionp.
const std::shared_ptr< LocalRegions::ExpansionVector >	GetExp () const
	This function returns the vector of elements in the expansion.
LocalRegions::ExpansionSharedPtr &	GetExp (int n) const
	This function returns (a shared pointer to) the local elemental expansion of the \(n^{\mathrm{th}}\) element.
LocalRegions::ExpansionSharedPtr &	GetExpFromGeomId (int n)
	This function returns (a shared pointer to) the local elemental expansion of the \(n^{\mathrm{th}}\) element given a global geometry ID.
LocalRegions::ExpansionSharedPtr &	GetExp (const Array< OneD, const NekDouble > &gloCoord)
	This function returns (a shared pointer to) the local elemental expansion containing the arbitrary point given by gloCoord.
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.
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.
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.
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.
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.
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.
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.
void	SetBndCondBwdWeight (const int index, const NekDouble value)
	Set the weight value for boundary conditions.
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.
void	GetBwdWeight (Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
	Get the weight value for boundary conditions for boundary average and jump calculations.
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.
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.
void	PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
	Copy and fill the Periodic boundaries.
void	PeriodicBwdRot (Array< OneD, Array< OneD, NekDouble > > &Bwd)
	Rotate Bwd trace for rotational periodicity boundaries when the flow is perpendicular to the rotation axis.
void	PeriodicDeriveBwdRot (TensorOfArray3D< NekDouble > &Bwd)
	Rotate Bwd trace derivative for rotational periodicity boundaries when the flow is perpendicular to the rotation axis.
void	RotLocalBwdTrace (Array< OneD, Array< OneD, NekDouble > > &Bwd)
	Rotate local Bwd trace across a rotational interface when the flow is perpendicular to the rotation axis.
void	RotLocalBwdDeriveTrace (TensorOfArray3D< NekDouble > &Bwd)
	Rotate local Bwd trace derivatives across a rotational interface when the flow is perpendicular to the rotation axis.
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, bool gridVelocity=false)
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	SetBCsToHomogeneous (void)
	Set boundary conditions to be homogeneous.
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.
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, RobinBCInfoSharedPtr >	GetRobinBCInfo ()
void	GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces=NullPeriodicMap)
std::vector< LibUtilities::FieldDefinitionsSharedPtr >	GetFieldDefinitions ()
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.
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.
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.
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.
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 >())
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)
std::shared_ptr< ExpList >	GetSharedThisPtr ()
	Returns a shared pointer to the current object.
std::shared_ptr< LibUtilities::SessionReader >	GetSession () const
	Returns the session object.
std::shared_ptr< LibUtilities::Comm >	GetComm () const
	Returns the comm object.
SpatialDomains::MeshGraphSharedPtr	GetGraph ()
LibUtilities::BasisSharedPtr	GetHomogeneousBasis (void)
std::shared_ptr< ExpList > &	GetPlane (int n)
const std::shared_ptr< GJPStabilisation >	GetGJPData (void)
void	CreateCollections (Collections::ImplementationType ImpType=Collections::eNoImpType)
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.
void	AddTraceJacToElmtJac (const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
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 LocTraceToTraceMapSharedPtr &	GetLocTraceToTraceMap () 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.
int	GetElmtToExpId (int elmtId)
	This function returns the index inside m_exp for a given geom id.
void	MultiplyByBlockMatrix (const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Transpose=false)
const Collections::CollectionVector &	GetCollections () const
	This function returns collections.
int	Get_coll_coeff_offset (int n) const
int	Get_coll_phys_offset (int n) const
const Array< OneD, const Array< OneD, NekDouble > > &	GetGridVelocity ()

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, const Collections::ImplementationType ImpType=Collections::eNoImpType)
	Discretises the boundary conditions.
void	GenerateBoundaryConditionExpansion (const Array< OneD, const MultiRegions::ExpListSharedPtr > &In, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
	Make copy of boundary conditions.
void	FindPeriodicTraces (const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
	Generate a associative map of periodic vertices in a mesh.
void	SetUpDG (const std::string="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
	Set up all DG member variables and maps.
ExpListSharedPtr &	v_GetTrace () override
AssemblyMapDGSharedPtr &	v_GetTraceMap (void) override
InterfaceMapDGSharedPtr &	v_GetInterfaceMap (void) override
const LocTraceToTraceMapSharedPtr &	v_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.
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.
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, bool gridVelocity=false) override
	This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray.
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, RobinBCInfoSharedPtr >	v_GetRobinBCInfo () override
const Array< OneD, const MultiRegions::ExpListSharedPtr > &	v_GetBndCondExpansions () override
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &	v_GetBndConditions () override
MultiRegions::ExpListSharedPtr &	v_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.
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..
void	v_SetBCsToHomogeneous (void) override
	Set boundary conditions to be homogeneous.
GlobalLinSysKey	v_HelmSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const StdRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing) override
	Solve the Helmholtz equation.
void	v_PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void	v_PeriodicBwdRot (Array< OneD, Array< OneD, NekDouble > > &Bwd) override
void	v_PeriodicDeriveBwdRot (TensorOfArray3D< NekDouble > &Bwd) override
void	v_RotLocalBwdTrace (Array< OneD, Array< OneD, NekDouble > > &Bwd) override
void	v_RotLocalBwdDeriveTrace (TensorOfArray3D< NekDouble > &Bwd) override
void	v_FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp) override
	Fill the weight with m_bndCondBndWeight.
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.
void	v_AddTraceIntegralToOffDiag (const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray) override
void	Rotate (Array< OneD, Array< OneD, NekDouble > > &Bwd, const int dir, const NekDouble angle, const int offset, const int npts)
	Rotate the slice [offset, offset + npts) of the 3D vector field in Bwd around the axis 'dir' by 'angle'.
void	DeriveRotate (TensorOfArray3D< NekDouble > &Bwd, const int dir, const NekDouble angle, const int offset, const int npts)
	Rotate the slice [offset, offset + npts) of the 3D vector field in Bwd around the axis 'dir' by 'angle'.
Protected Member Functions inherited from Nektar::MultiRegions::ExpList
std::shared_ptr< DNekMat >	GenGlobalMatrixFull (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.
const DNekScalBlkMatSharedPtr &	GetBlockMatrix (const GlobalMatrixKey &gkey)
std::shared_ptr< GlobalMatrix >	GenGlobalMatrix (const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
	Generates a global matrix from the given key and map.
void	GlobalEigenSystem (const std::shared_ptr< DNekMat > &Gmat, Array< OneD, NekDouble > &EigValsReal, Array< OneD, NekDouble > &EigValsImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
std::shared_ptr< GlobalLinSys >	GenGlobalLinSys (const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
	This operation constructs the global linear system of type mkey.
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 LocToGloBaseMap.
virtual size_t	v_GetNumElmts (void)
virtual void	v_Upwind (const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
virtual void	v_Upwind (const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
virtual const Array< OneD, const int > &	v_GetTraceBndMap ()
virtual void	v_GetNormals (Array< OneD, Array< OneD, NekDouble > > &normals)
	Populate normals with the normals of all expansions.
virtual void	v_AddTraceQuadPhysToOffDiag (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
virtual const std::vector< bool > &	v_GetLeftAdjacentFaces (void) const
virtual void	v_MultiplyByInvMassMatrix (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual GlobalLinSysKey	v_LinearAdvectionDiffusionReactionSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const StdRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
virtual GlobalLinSysKey	v_LinearAdvectionReactionSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const StdRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
virtual void	v_ImposeDirichletConditions (Array< OneD, NekDouble > &outarray)
virtual void	v_ImposeNeumannConditions (Array< OneD, NekDouble > &outarray)
virtual void	v_ImposeRobinConditions (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_AvgAssemble (bool UseComm)
virtual void	v_AvgAssemble (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseComm)
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 ()
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::FieldDefinitionsSharedPtr >	v_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)
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 NekDouble >	v_HomogeneousEnergy (void)
virtual LibUtilities::TranspositionSharedPtr	v_GetTransposition (void)
virtual NekDouble	v_GetHomoLen (void)
virtual void	v_SetHomoLen (const NekDouble lhom)
virtual Array< OneD, const unsigned int >	v_GetZIDs (void)
virtual Array< OneD, const unsigned int >	v_GetYIDs (void)
virtual void	v_PhysInterp1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void	v_PhysGalerkinProjection1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void	v_ClearGlobalLinSysManager (void)
virtual int	v_GetPoolCount (std::string)
virtual void	v_UnsetGlobalLinSys (GlobalLinSysKey, bool)
virtual LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > &	v_GetGlobalLinSysManager (void)
void	ExtractFileBCs (const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
virtual LibUtilities::BasisSharedPtr	v_GetHomogeneousBasis (void)
virtual void	v_SetHomo1DSpecVanVisc (Array< OneD, NekDouble > visc)
virtual std::shared_ptr< ExpList > &	v_GetPlane (int n)
virtual const std::shared_ptr< GJPStabilisation >	v_GetGJPData (void)

Protected Attributes
Array< OneD, SpatialDomains::BoundaryConditionShPtr >	m_bndConditions
	An array which contains the information about the boundary condition structure definition on the different boundary regions.
Array< OneD, MultiRegions::ExpListSharedPtr >	m_bndCondExpansions
	An object which contains the discretised boundary conditions.
Array< OneD, NekDouble >	m_bndCondBndWeight
InterfaceMapDGSharedPtr	m_interfaceMap
	Interfaces mapping for trace space.
GlobalLinSysMapShPtr	m_globalBndMat
	Global boundary matrix.
ExpListSharedPtr	m_trace
	Trace space storage for points between elements.
MultiRegions::ExpListSharedPtr	m_locElmtTrace
	Local Elemental trace expansions.
AssemblyMapDGSharedPtr	m_traceMap
	Local to global DG mapping for trace space.
std::set< int >	m_boundaryTraces
	A set storing the global IDs of any boundary Verts.
PeriodicMap	m_periodicVerts
	A map which identifies groups of periodic vertices.
PeriodicMap	m_periodicEdges
	A map which identifies pairs of periodic edges.
PeriodicMap	m_periodicFaces
	A map which identifies pairs of periodic faces.
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.
std::vector< int >	m_periodicBwdCopy
std::vector< bool >	m_leftAdjacentTraces
LocTraceToTraceMapSharedPtr	m_locTraceToTraceMap
Protected Attributes inherited from Nektar::MultiRegions::ExpList
SpatialDomains::EntityHolder1D	m_holder
	Pointer holder for PulseWaveSolver.
ExpansionType	m_expType
	Expansion type.
LibUtilities::CommSharedPtr	m_comm
	Communicator.
LibUtilities::SessionReaderSharedPtr	m_session
	Session.
SpatialDomains::MeshGraphSharedPtr	m_graph
	Mesh associated with this expansion list.
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\).
int	m_npoints
Array< OneD, NekDouble >	m_coeffs
	Concatenation of all local expansion coefficients.
Array< OneD, NekDouble >	m_phys
	The global expansion evaluated at the quadrature points.
bool	m_physState
	The state of the array m_phys.
std::shared_ptr< LocalRegions::ExpansionVector >	m_exp
	The list of local expansions.
Collections::CollectionVector	m_collections
std::vector< bool >	m_collectionsDoInit
	Vector of bools to act as an initialise on first call flag.
std::vector< int >	m_coll_coeff_offset
	Offset of elemental data into the array m_coeffs.
std::vector< int >	m_coll_phys_offset
	Offset of elemental data into the array m_phys.
Array< OneD, int >	m_coeff_offset
	Offset of elemental data into the array m_coeffs.
Array< OneD, int >	m_phys_offset
	Offset of elemental data into the array m_phys.
Array< OneD, std::pair< int, int > >	m_coeffsToElmt
	m_coeffs to elemental value map
BlockMatrixMapShPtr	m_blockMat
bool	m_WaveSpace
std::unordered_map< int, int >	m_elmtToExpId
	Mapping from geometry ID of element to index inside m_exp.
Array< OneD, Array< OneD, NekDouble > >	m_gridVelocity
	Grid velocity at quadrature points.

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

    : ExpList(), m_bndConditions(), m_bndCondExpansions(),
      m_trace(NullExpListSharedPtr)
{
}

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

graph1D	A mesh, containing information about the domain and the spectral/hp element expansions.
bcs	Information about the enforced boundary conditions.
variable	The session variable associated with the boundary conditions to enforce.
solnType	Type of global system to use.

Definition at line 81 of file DisContField.cpp.

    : ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
      m_trace(NullExpListSharedPtr)
{
    std::string bcvar;
    if (bcvariable == "NotSet")
    {
        bcvar = variable;
    }
    else
    {
        bcvar = bcvariable;
    }
 
    if (bcvar.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditions bcs(m_session, graph);
 
        GenerateBoundaryConditionExpansion(graph, bcs, bcvar,
                                           DeclareCoeffPhysArrays, ImpType);
        if (DeclareCoeffPhysArrays)
        {
            EvaluateBoundaryConditions(0.0, bcvar);
        }
        ApplyGeomInfo();
        FindPeriodicTraces(bcs, bcvar);
    }
 
    int i, cnt;
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        MultiRegions::ExpListSharedPtr locExpList;
        locExpList = m_bndCondExpansions[i];
 
        for (int e = 0; e < locExpList->GetExpSize(); ++e)
        {
            (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
                                                   locExpList->GetExp(e));
            locExpList->GetExp(e)->SetAdjacentElementExp(
                TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
        }
        cnt += m_bndCondExpansions[i]->GetExpSize();
    }
 
    if (SetUpJustDG)
    {
        SetUpDG(variable, ImpType);
    }
    else
    {
        if (m_session->DefinesSolverInfo("PROJECTION"))
        {
            std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
            if ((ProjectStr == "MixedCGDG") ||
                (ProjectStr == "Mixed_CG_Discontinuous"))
            {
                SetUpDG(variable);
            }
            else
            {
                SetUpPhysNormals();
            }
        }
        else
        {
            SetUpPhysNormals();
        }
    }
}

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

domain

Subdomain specified in the inputfile from which the DisContField is set up

Definition at line 608 of file DisContField.cpp.

    : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
              comm, ImpType)
{
    if (variable.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditionsSharedPtr DomBCs =
            GetDomainBCs(domain, Allbcs, variable);
 
        GenerateBoundaryConditionExpansion(m_graph, *DomBCs, variable, true,
                                           ImpType);
        EvaluateBoundaryConditions(0.0, variable);
        ApplyGeomInfo();
        FindPeriodicTraces(*DomBCs, variable);
    }
 
    SetUpDG(variable);
}

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

In	Existing DisContField object to copy.

Definition at line 638 of file DisContField.cpp.

    : ExpList(In, DeclareCoeffPhysArrays), m_bndConditions(In.m_bndConditions),
      m_bndCondExpansions(In.m_bndCondExpansions),
      m_bndCondBndWeight(In.m_bndCondBndWeight),
      m_interfaceMap(In.m_interfaceMap), m_globalBndMat(In.m_globalBndMat),
      m_traceMap(In.m_traceMap), m_boundaryTraces(In.m_boundaryTraces),
      m_periodicVerts(In.m_periodicVerts),
      m_periodicFwdCopy(In.m_periodicFwdCopy),
      m_periodicBwdCopy(In.m_periodicBwdCopy),
      m_leftAdjacentTraces(In.m_leftAdjacentTraces),
      m_locTraceToTraceMap(In.m_locTraceToTraceMap)
{
    if (In.m_trace) // independent trace space
    {
        m_trace = MemoryManager<ExpList>::AllocateSharedPtr(
            *In.m_trace, DeclareCoeffPhysArrays);
    }
 
    if (In.m_locElmtTrace)
    {
        m_locElmtTrace = MemoryManager<ExpList>::AllocateSharedPtr(
            *In.m_locElmtTrace, DeclareCoeffPhysArrays);
    }
}

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

DisContField() [5/6]

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

Definition at line 668 of file DisContField.cpp.

    : ExpList(In, DeclareCoeffPhysArrays)
{
    m_trace = NullExpListSharedPtr;
 
    // Set up boundary conditions for this variable.
    // Do not set up BCs if default variable
    if (variable.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditions bcs(m_session, graph);
        GenerateBoundaryConditionExpansion(In.m_bndCondExpansions, bcs,
                                           variable);
 
        if (DeclareCoeffPhysArrays)
        {
            EvaluateBoundaryConditions(0.0, variable);
        }
 
        if (!SameTypeOfBoundaryConditions(In))
        {
            // Find periodic edges for this variable.
            FindPeriodicTraces(bcs, variable);
 
            if (SetUpJustDG)
            {
                SetUpDG(variable);
            }
            else
            {
                // set elmt edges to point to robin bc edges if required
                int i, cnt = 0;
                Array<OneD, int> ElmtID, TraceID;
                GetBoundaryToElmtMap(ElmtID, TraceID);
 
                for (i = 0; i < m_bndCondExpansions.size(); ++i)
                {
                    MultiRegions::ExpListSharedPtr locExpList;
 
                    int e;
                    locExpList = m_bndCondExpansions[i];
 
                    for (e = 0; e < locExpList->GetExpSize(); ++e)
                    {
                        (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
                            TraceID[cnt + e], locExpList->GetExp(e));
                        locExpList->GetExp(e)->SetAdjacentElementExp(
                            TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
                    }
 
                    cnt += m_bndCondExpansions[i]->GetExpSize();
                }
 
                if (m_session->DefinesSolverInfo("PROJECTION"))
                {
                    std::string ProjectStr =
                        m_session->GetSolverInfo("PROJECTION");
 
                    if ((ProjectStr == "MixedCGDG") ||
                        (ProjectStr == "Mixed_CG_Discontinuous"))
                    {
                        SetUpDG(variable);
                    }
                    else
                    {
                        SetUpPhysNormals();
                    }
                }
                else
                {
                    SetUpPhysNormals();
                }
            }
        }
        else
        {
            m_globalBndMat       = In.m_globalBndMat;
            m_trace              = In.m_trace;
            m_locElmtTrace       = In.m_locElmtTrace;
            m_traceMap           = In.m_traceMap;
            m_interfaceMap       = In.m_interfaceMap;
            m_locTraceToTraceMap = In.m_locTraceToTraceMap;
            m_periodicVerts      = In.m_periodicVerts;
            m_periodicEdges      = In.m_periodicEdges;
            m_periodicFaces      = In.m_periodicFaces;
            m_periodicFwdCopy    = In.m_periodicFwdCopy;
            m_periodicBwdCopy    = In.m_periodicBwdCopy;
            m_boundaryTraces     = In.m_boundaryTraces;
            m_leftAdjacentTraces = In.m_leftAdjacentTraces;
 
            if (!SetUpJustDG)
            {
                // set elmt edges to point to robin bc edges if required
                int i, cnt = 0;
                Array<OneD, int> ElmtID, TraceID;
                GetBoundaryToElmtMap(ElmtID, TraceID);
 
                for (i = 0; i < m_bndCondExpansions.size(); ++i)
                {
                    MultiRegions::ExpListSharedPtr locExpList;
 
                    int e;
                    locExpList = m_bndCondExpansions[i];
 
                    for (e = 0; e < locExpList->GetExpSize(); ++e)
                    {
                        (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
                            TraceID[cnt + e], locExpList->GetExp(e));
                        locExpList->GetExp(e)->SetAdjacentElementExp(
                            TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
                    }
                    cnt += m_bndCondExpansions[i]->GetExpSize();
                }
 
                SetUpPhysNormals();
            }
        }
    }
}

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_locElmtTrace, 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

In	Existing ExpList object to copy.

Definition at line 794 of file DisContField.cpp.

                                            : ExpList(In)
{
}

~DisContField()

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

override

Destructor.

Definition at line 801 of file DisContField.cpp.

{
}

Member Function Documentation

DeriveRotate()

void Nektar::MultiRegions::DisContField::DeriveRotate ( TensorOfArray3D< NekDouble > &  Bwd, const int  dir, const NekDouble  angle, const int  offset, const int  npts  )

protected

Rotate the slice [offset, offset + npts) of the 3D vector field in Bwd around the axis 'dir' by 'angle'.

Definition at line 4962 of file DisContField.cpp.

{
    // Precompute trigonometric values using correct angle input
    const NekDouble cosA = cos(angle);
    const NekDouble sinA = sin(angle);
 
    // Define rotation matrix R (3x3) using Array<OneD, NekDouble> in
    // **column-major order**
    Array<OneD, NekDouble> R(9, 0.0);
 
    if (dir == 0) // Rotation around X-axis
    {
        R[0] = 1.0;
        R[3] = 0.0;
        R[6] = 0.0;
        R[1] = 0.0;
        R[4] = cosA;
        R[7] = -sinA;
        R[2] = 0.0;
        R[5] = sinA;
        R[8] = cosA;
    }
    else if (dir == 1) // Rotation around Y-axis
    {
        R[0] = cosA;
        R[3] = 0.0;
        R[6] = sinA;
        R[1] = 0.0;
        R[4] = 1.0;
        R[7] = 0.0;
        R[2] = -sinA;
        R[5] = 0.0;
        R[8] = cosA;
    }
    else if (dir == 2) // Rotation around Z-axis
    {
        R[0] = cosA;
        R[3] = -sinA;
        R[6] = 0.0;
        R[1] = sinA;
        R[4] = cosA;
        R[7] = 0.0;
        R[2] = 0.0;
        R[5] = 0.0;
        R[8] = 1.0;
    }
    else
    {
        NEKERROR(
            ErrorUtil::efatal,
            "Invalid rotation axis for first-order derivative transformation.");
    }
 
    for (int p = offset; p < offset + npts; ++p)
    {
        // Allocate gradient tensor (3x3) in **column-major order**
        Array<OneD, NekDouble> gradU(9, 0.0);
 
        // Load original velocity gradient tensor into **column-major** format
        gradU[0] = Bwd[0][1][p]; // ∂u_x/∂x
        gradU[3] = Bwd[1][1][p]; // ∂u_x/∂y
        gradU[6] = Bwd[2][1][p]; // ∂u_x/∂z
 
        gradU[1] = Bwd[0][2][p]; // ∂u_y/∂x
        gradU[4] = Bwd[1][2][p]; // ∂u_y/∂y
        gradU[7] = Bwd[2][2][p]; // ∂u_y/∂z
 
        gradU[2] = Bwd[0][3][p]; // ∂u_z/∂x
        gradU[5] = Bwd[1][3][p]; // ∂u_z/∂y
        gradU[8] = Bwd[2][3][p]; // ∂u_z/∂z
 
        // Allocate intermediate matrix R * gradU in **column-major order**
        Array<OneD, NekDouble> R_gradU(9, 0.0);
 
        // Perform matrix multiplication R * gradU using Blas::Dgemm in
        // **column-major order**
        Blas::Dgemm('N', 'N', 3, 3, 3, 1.0, R.data(), 3, gradU.data(), 3, 0.0,
                    R_gradU.data(), 3);
 
        // Allocate final rotated gradient (R * gradU) * R^T in **column-major
        // order**
        Array<OneD, NekDouble> gradU_rotated(9, 0.0);
 
        // Perform matrix multiplication (R * gradU) * R^T using Blas::Dgemm in
        // **column-major order**
        Blas::Dgemm('N', 'T', 3, 3, 3, 1.0, R_gradU.data(), 3, R.data(), 3, 0.0,
                    gradU_rotated.data(), 3);
 
        // Store rotated gradients back in Bwd in **column-major order**
        Bwd[0][1][p] = gradU_rotated[0]; // ∂u_x'/∂x'
        Bwd[1][1][p] = gradU_rotated[3]; // ∂u_x'/∂y'
        Bwd[2][1][p] = gradU_rotated[6]; // ∂u_x'/∂z'
 
        Bwd[0][2][p] = gradU_rotated[1]; // ∂u_y'/∂x'
        Bwd[1][2][p] = gradU_rotated[4]; // ∂u_y'/∂y'
        Bwd[2][2][p] = gradU_rotated[7]; // ∂u_y'/∂z'
 
        Bwd[0][3][p] = gradU_rotated[2]; // ∂u_z'/∂x'
        Bwd[1][3][p] = gradU_rotated[5]; // ∂u_z'/∂y'
        Bwd[2][3][p] = gradU_rotated[8]; // ∂u_z'/∂z'
    }
}

References Blas::Dgemm(), Nektar::ErrorUtil::efatal, and NEKERROR.

Referenced by v_PeriodicDeriveBwdRot().

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

outarray

The resulting field \( u^* \).

Definition at line 4652 of file DisContField.cpp.

{
    int i, cnt, e, ncoeff_trace;
    Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    StdRegions::Orientation edgedir;
 
    int nq_elmt, nm_elmt;
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
    Array<OneD, NekDouble> tmp_coeffs;
    m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
 
    trace_lambda = loc_lambda;
 
    int dim = (m_expType == e2D) ? 2 : 3;
 
    int num_points[] = {0, 0, 0};
    int num_modes[]  = {0, 0, 0};
 
    // Calculate Q using standard DG formulation.
    for (i = cnt = 0; i < GetExpSize(); ++i)
    {
        nq_elmt = (*m_exp)[i]->GetTotPoints();
        nm_elmt = (*m_exp)[i]->GetNcoeffs();
        qrhs    = Array<OneD, NekDouble>(nq_elmt);
        qrhs1   = Array<OneD, NekDouble>(nq_elmt);
        force   = Array<OneD, NekDouble>(2 * nm_elmt);
        out_tmp = force + nm_elmt;
        LocalRegions::ExpansionSharedPtr ppExp;
 
        for (int j = 0; j < dim; ++j)
        {
            num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
            num_modes[j]  = (*m_exp)[i]->GetBasis(j)->GetNumModes();
        }
 
        // Probably a better way of setting up lambda than this.  Note
        // cannot use PutCoeffsInToElmts since lambda space is mapped
        // during the solve.
        int nTraces = (*m_exp)[i]->GetNtraces();
        Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
 
        for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
        {
            edgedir        = (*m_exp)[i]->GetTraceOrient(e);
            ncoeff_trace   = elmtToTrace[i][e]->GetNcoeffs();
            traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
            Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
            if (dim == 2)
            {
                elmtToTrace[i][e]->SetCoeffsToOrientation(
                    edgedir, traceCoeffs[e], traceCoeffs[e]);
            }
            else
            {
                (*m_exp)[i]
                    ->as<LocalRegions::Expansion3D>()
                    ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
            }
            trace_lambda = trace_lambda + ncoeff_trace;
        }
 
        // creating orthogonal expansion (checking if we have quads or
        // triangles)
        LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
        switch (shape)
        {
            case LibUtilities::eQuadrilateral:
            {
                const LibUtilities::PointsKey PkeyQ1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyQ2(
                    num_points[1], LibUtilities::eGaussLobattoLegendre);
                LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyQ1);
                LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
                                              num_modes[1], PkeyQ2);
                SpatialDomains::QuadGeom *qGeom =
                    dynamic_cast<SpatialDomains::QuadGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
                    BkeyQ1, BkeyQ2, qGeom);
            }
            break;
            case LibUtilities::eTriangle:
            {
                const LibUtilities::PointsKey PkeyT1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyT2(
                    num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
                LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyT1);
                LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
                                              num_modes[1], PkeyT2);
                SpatialDomains::TriGeom *tGeom =
                    dynamic_cast<SpatialDomains::TriGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(
                    BkeyT1, BkeyT2, tGeom);
            }
            break;
            case LibUtilities::eHexahedron:
            {
                const LibUtilities::PointsKey PkeyH1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyH2(
                    num_points[1], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyH3(
                    num_points[2], LibUtilities::eGaussLobattoLegendre);
                LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyH1);
                LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
                                              num_modes[1], PkeyH2);
                LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
                                              num_modes[2], PkeyH3);
                SpatialDomains::HexGeom *hGeom =
                    dynamic_cast<SpatialDomains::HexGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
                    BkeyH1, BkeyH2, BkeyH3, hGeom);
            }
            break;
            case LibUtilities::eTetrahedron:
            {
                const LibUtilities::PointsKey PkeyT1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyT2(
                    num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
                const LibUtilities::PointsKey PkeyT3(
                    num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
                LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyT1);
                LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
                                              num_modes[1], PkeyT2);
                LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
                                              num_modes[2], PkeyT3);
                SpatialDomains::TetGeom *tGeom =
                    dynamic_cast<SpatialDomains::TetGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
                    BkeyT1, BkeyT2, BkeyT3, tGeom);
            }
            break;
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Wrong shape type, HDG postprocessing is not "
                         "implemented");
        };
 
        // DGDeriv
        // (d/dx w, d/dx q_0)
        (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], traceCoeffs, out_tmp);
        (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
        ppExp->IProductWRTDerivBase(0, qrhs, force);
 
        // + (d/dy w, d/dy q_1)
        (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], traceCoeffs, out_tmp);
 
        (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
        ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
 
        Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
 
        // determine force[0] = (1,u)
        (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
        force[0] = (*m_exp)[i]->Integral(qrhs);
 
        // multiply by inverse Laplacian matrix
        // get matrix inverse
        LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
                                       ppExp->DetShapeType(), *ppExp);
        DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
 
        NekVector<NekDouble> in(nm_elmt, force, eWrapper);
        NekVector<NekDouble> out(nm_elmt);
 
        out = (*lapsys) * in;
 
        // Transforming back to modified basis
        Array<OneD, NekDouble> work(nq_elmt);
        ppExp->BwdTrans(out.GetPtr(), work);
        (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
    }
}

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

{
 
    // Fill boundary conditions into missing elements
    if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
    {
        // Fill boundary conditions into missing elements
        for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
        {
            if (bndConditions[n]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
 
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann ||
                     bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() !=
                     SpatialDomains::ePeriodic)
            {
                ASSERTL0(false, "Method not set up for this "
                                "boundary condition.");
            }
        }
    }
    else
    {
        for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
        {
            if (bndConditions[n]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id1  = bndCondExpansions[n]->GetPhys_Offset(e);
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
 
                    Vmath::Vcopy(npts, &(bndCondExpansions[n]->GetPhys())[id1],
                                 1, &Bwd[id2], 1);
                }
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann ||
                     bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
            {
                auto ne = m_bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id1  = bndCondExpansions[n]->GetPhys_Offset(e);
 
                    ASSERTL0((bndCondExpansions[n]->GetPhys())[id1] == 0.0,
                             "Method not set up for non-zero "
                             "Neumann boundary condition");
                    auto id2 = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
 
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
 
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eNotDefined)
            {
                // Do nothing
            }
            else if (bndConditions[n]->GetBoundaryConditionType() !=
                     SpatialDomains::ePeriodic)
            {
                NEKERROR(ErrorUtil::efatal,
                         "Method not set up for this boundary "
                         "condition.");
            }
        }
    }
}

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

bcs	Information about the boundary conditions.
variable	Specifies the field.

Definition at line 1003 of file DisContField.cpp.

{
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
 
    switch (m_expType)
    {
        case e1D:
        {
            int i, region1ID, region2ID;
 
            SpatialDomains::BoundaryConditionShPtr locBCond;
            map<int, int> BregionToVertMap;
 
            // Construct list of all periodic Region and their
            // global vertex on this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
                int id =
                    it.second->begin()->second->m_geomVec[0]->GetGlobalID();
 
                BregionToVertMap[it.first] = id;
            }
 
            set<int> islocal;
 
            int n = vComm->GetSize();
            int p = vComm->GetRank();
 
            Array<OneD, int> nregions(n, 0);
            nregions[p] = BregionToVertMap.size();
            vComm->AllReduce(nregions, LibUtilities::ReduceSum);
 
            int totRegions = Vmath::Vsum(n, nregions, 1);
 
            Array<OneD, int> regOffset(n, 0);
 
            for (i = 1; i < n; ++i)
            {
                regOffset[i] = regOffset[i - 1] + nregions[i - 1];
            }
 
            Array<OneD, int> bregmap(totRegions, 0);
            Array<OneD, int> bregid(totRegions, 0);
 
            i = regOffset[p];
            for (auto &iit : BregionToVertMap)
            {
                bregid[i]    = iit.first;
                bregmap[i++] = iit.second;
                islocal.insert(iit.first);
            }
 
            vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
            vComm->AllReduce(bregid, LibUtilities::ReduceSum);
 
            for (int i = 0; i < totRegions; ++i)
            {
                BregionToVertMap[bregid[i]] = bregmap[i];
            }
 
            // Construct list of all periodic pairs local to this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                ASSERTL0(BregionToVertMap.count(region1ID) != 0,
                         "Cannot determine vertex of region1ID");
 
                ASSERTL0(BregionToVertMap.count(region2ID) != 0,
                         "Cannot determine vertex of region2ID");
 
                PeriodicEntity ent(BregionToVertMap[region2ID],
                                   StdRegions::eNoOrientation,
                                   islocal.count(region2ID) != 0);
 
                m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
            }
        }
        break;
        case e2D:
        {
            int region1ID, region2ID, i, j, k, cnt;
            SpatialDomains::BoundaryConditionShPtr locBCond;
 
            SpatialDomains::CompositeOrdering compOrder =
                m_graph->GetCompositeOrdering();
            SpatialDomains::BndRegionOrdering bndRegOrder =
                m_graph->GetBndRegionOrdering();
            SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
 
            // Unique collection of pairs of periodic composites
            // (i.e. if composites 1 and 2 are periodic then this
            // map will contain either the pair (1,2) or (2,1) but
            // not both).
            map<int, int> perComps;
            map<int, vector<int>> allVerts;
            set<int> locVerts;
            map<int, pair<int, StdRegions::Orientation>> allEdges;
 
            // Set up a set of all local verts and edges.
            for (i = 0; i < (*m_exp).size(); ++i)
            {
                for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetVid(j);
                    locVerts.insert(id);
                }
            }
 
            // Construct list of all periodic pairs local to this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                // From this identify composites. Note that in
                // serial this will be an empty map.
                int cId1, cId2;
                if (vComm->GetSize() == 1)
                {
                    cId1 = it.second->begin()->first;
                    cId2 = bregions.find(region2ID)->second->begin()->first;
                }
                else
                {
                    cId1 = bndRegOrder.find(region1ID)->second[0];
                    cId2 = bndRegOrder.find(region2ID)->second[0];
                }
 
                ASSERTL0(it.second->size() == 1,
                         "Boundary region " + std::to_string(region1ID) +
                             " should only contain 1 composite.");
 
                vector<unsigned int> tmpOrder;
 
                // Construct set containing all periodic edgesd on
                // this process
                SpatialDomains::CompositeSharedPtr c =
                    it.second->begin()->second;
 
                for (i = 0; i < c->m_geomVec.size(); ++i)
                {
                    SpatialDomains::SegGeom *segGeom =
                        dynamic_cast<SpatialDomains::SegGeom *>(
                            c->m_geomVec[i]);
                    ASSERTL0(segGeom, "Unable to cast to SegGeom*");
 
                    SpatialDomains::GeometryLinkSharedPtr elmt =
                        m_graph->GetElementsFromEdge(segGeom);
                    ASSERTL0(elmt->size() == 1,
                             "The periodic boundaries belong to "
                             "more than one element of the mesh");
 
                    SpatialDomains::Geometry2D *geom =
                        dynamic_cast<SpatialDomains::Geometry2D *>(
                            elmt->at(0).first);
 
                    allEdges[c->m_geomVec[i]->GetGlobalID()] =
                        make_pair(elmt->at(0).second,
                                  geom->GetEorient(elmt->at(0).second));
 
                    // In serial mesh partitioning will not have occurred so
                    // need to fill composite ordering map manually.
                    if (vComm->GetSize() == 1)
                    {
                        tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
                    }
 
                    vector<int> vertList(2);
                    vertList[0] = segGeom->GetVid(0);
                    vertList[1] = segGeom->GetVid(1);
                    allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
                }
 
                if (vComm->GetSize() == 1)
                {
                    compOrder[it.second->begin()->first] = tmpOrder;
                }
 
                // See if we already have either region1 or
                // region2 stored in perComps map.
                if (perComps.count(cId1) == 0)
                {
                    if (perComps.count(cId2) == 0)
                    {
                        perComps[cId1] = cId2;
                    }
                    else
                    {
                        std::stringstream ss;
                        ss << "Boundary region " << cId2 << " should be "
                           << "periodic with " << perComps[cId2] << " but "
                           << "found " << cId1 << " instead!";
                        ASSERTL0(perComps[cId2] == cId1, ss.str());
                    }
                }
                else
                {
                    std::stringstream ss;
                    ss << "Boundary region " << cId1 << " should be "
                       << "periodic with " << perComps[cId1] << " but "
                       << "found " << cId2 << " instead!";
                    ASSERTL0(perComps[cId1] == cId1, ss.str());
                }
            }
 
            // In case of periodic partition being split over many composites
            // may not have all periodic matches so check this
            int idmax = -1;
            for (auto &cIt : perComps)
            {
                idmax = max(idmax, cIt.first);
                idmax = max(idmax, cIt.second);
            }
            vComm->AllReduce(idmax, LibUtilities::ReduceMax);
            idmax++;
            Array<OneD, int> perid(idmax, -1);
            for (auto &cIt : perComps)
            {
                perid[cIt.first] = cIt.second;
            }
            vComm->AllReduce(perid, LibUtilities::ReduceMax);
            // update all partitions
            for (int i = 0; i < idmax; ++i)
            {
                if (perid[i] > -1)
                {
                    // skip if complemetary relationship has
                    // already been speficied in list
                    if (perComps.count(perid[i]))
                    {
                        continue;
                    }
                    perComps[i] = perid[i];
                }
            }
 
            // Process local edge list to obtain relative edge
            // orientations.
            int n = vComm->GetSize();
            int p = vComm->GetRank();
            int totEdges;
            Array<OneD, int> edgecounts(n, 0);
            Array<OneD, int> edgeoffset(n, 0);
            Array<OneD, int> vertoffset(n, 0);
 
            edgecounts[p] = allEdges.size();
            vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
 
            edgeoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
            }
 
            totEdges = Vmath::Vsum(n, edgecounts, 1);
            Array<OneD, int> edgeIds(totEdges, 0);
            Array<OneD, int> edgeIdx(totEdges, 0);
            Array<OneD, int> edgeOrient(totEdges, 0);
            Array<OneD, int> edgeVerts(totEdges, 0);
 
            auto sIt = allEdges.begin();
 
            for (i = 0; sIt != allEdges.end(); ++sIt)
            {
                edgeIds[edgeoffset[p] + i]     = sIt->first;
                edgeIdx[edgeoffset[p] + i]     = sIt->second.first;
                edgeOrient[edgeoffset[p] + i]  = sIt->second.second;
                edgeVerts[edgeoffset[p] + i++] = allVerts[sIt->first].size();
            }
 
            vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeIdx, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
 
            // Calculate number of vertices on each processor.
            Array<OneD, int> procVerts(n, 0);
            int nTotVerts;
 
            // Note if there are no periodic edges at all calling Vsum will
            // cause a segfault.
            if (totEdges > 0)
            {
                nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
            }
            else
            {
                nTotVerts = 0;
            }
 
            for (i = 0; i < n; ++i)
            {
                if (edgecounts[i] > 0)
                {
                    procVerts[i] = Vmath::Vsum(edgecounts[i],
                                               edgeVerts + edgeoffset[i], 1);
                }
                else
                {
                    procVerts[i] = 0;
                }
            }
            vertoffset[0] = 0;
 
            for (i = 1; i < n; ++i)
            {
                vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
            }
 
            Array<OneD, int> traceIds(nTotVerts, 0);
            for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
            {
                for (j = 0; j < allVerts[sIt->first].size(); ++j)
                {
                    traceIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
                }
            }
 
            vComm->AllReduce(traceIds, LibUtilities::ReduceSum);
 
            // For simplicity's sake create a map of edge id -> orientation.
            map<int, StdRegions::Orientation> orientMap;
            map<int, vector<int>> vertMap;
 
            for (cnt = i = 0; i < totEdges; ++i)
            {
                ASSERTL0(orientMap.count(edgeIds[i]) == 0,
                         "Already found edge in orientation map!");
 
                // Work out relative orientations. To avoid having
                // to exchange vertex locations, we figure out if
                // the edges are backwards or forwards orientated
                // with respect to a forwards orientation that is
                // CCW. Since our local geometries are
                // forwards-orientated with respect to the
                // Cartesian axes, we need to invert the
                // orientation for the top and left edges of a
                // quad and the left edge of a triangle.
                StdRegions::Orientation o =
                    (StdRegions::Orientation)edgeOrient[i];
 
                if (edgeIdx[i] > 1)
                {
                    o = o == StdRegions::eForwards ? StdRegions::eBackwards
                                                   : StdRegions::eForwards;
                }
 
                orientMap[edgeIds[i]] = o;
 
                vector<int> verts(edgeVerts[i]);
 
                for (j = 0; j < edgeVerts[i]; ++j)
                {
                    verts[j] = traceIds[cnt++];
                }
                vertMap[edgeIds[i]] = verts;
            }
 
            // Go through list of composites and figure out which
            // edges are parallel from original ordering in
            // session file. This includes composites which are
            // not necessarily on this process.
            map<int, int> allCompPairs;
 
            // Store temporary map of periodic vertices which will hold all
            // periodic vertices on the entire mesh so that doubly periodic
            // vertices can be counted properly across partitions. Local
            // vertices are copied into m_periodicVerts at the end of the
            // function.
            PeriodicMap periodicVerts;
 
            for (auto &cIt : perComps)
            {
                SpatialDomains::CompositeSharedPtr c[2];
                const int id1    = cIt.first;
                const int id2    = cIt.second;
                std::string id1s = std::to_string(id1);
                std::string id2s = std::to_string(id2);
 
                if (compMap.count(id1) > 0)
                {
                    c[0] = compMap[id1];
                }
 
                if (compMap.count(id2) > 0)
                {
                    c[1] = compMap[id2];
                }
 
                // Loop over composite ordering to construct list of all
                // periodic edges regardless of whether they are on this
                // process.
                map<int, int> compPairs;
 
                ASSERTL0(compOrder.count(id1) > 0,
                         "Unable to find composite " + id1s + " in order map.");
                ASSERTL0(compOrder.count(id2) > 0,
                         "Unable to find composite " + id2s + " in order map.");
                WARNINGL0(
                    compOrder[id1].size() == compOrder[id2].size(),
                    "Periodic composites " + id1s + " and " + id2s +
                        " should have the same number of elements. Have " +
                        std::to_string(compOrder[id1].size()) + " vs " +
                        std::to_string(compOrder[id2].size()));
                ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
                                                        id1s + " and " + id2s +
                                                        " are empty!");
 
                // TODO: Add more checks.
                for (i = 0; i < compOrder[id1].size(); ++i)
                {
                    int eId1 = compOrder[id1][i];
                    int eId2 = compOrder[id2][i];
 
                    ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
 
                    if (compPairs.count(eId2) != 0)
                    {
                        ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
                    }
                    compPairs[eId1] = eId2;
                }
 
                // Construct set of all edges that we have
                // locally on this processor.
                set<int> locEdges;
                for (i = 0; i < 2; ++i)
                {
                    if (!c[i])
                    {
                        continue;
                    }
 
                    if (c[i]->m_geomVec.size() > 0)
                    {
                        for (j = 0; j < c[i]->m_geomVec.size(); ++j)
                        {
                            locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
                        }
                    }
                }
 
                // Loop over all edges in the geometry composite.
                for (auto &pIt : compPairs)
                {
                    int ids[2]    = {pIt.first, pIt.second};
                    bool local[2] = {locEdges.count(pIt.first) > 0,
                                     locEdges.count(pIt.second) > 0};
 
                    ASSERTL0(orientMap.count(ids[0]) > 0 &&
                                 orientMap.count(ids[1]) > 0,
                             "Unable to find edge in orientation map.");
 
                    allCompPairs[pIt.first]  = pIt.second;
                    allCompPairs[pIt.second] = pIt.first;
 
                    for (i = 0; i < 2; ++i)
                    {
                        if (!local[i])
                        {
                            continue;
                        }
 
                        int other = (i + 1) % 2;
 
                        StdRegions::Orientation o =
                            orientMap[ids[i]] == orientMap[ids[other]]
                                ? StdRegions::eBackwards
                                : StdRegions::eForwards;
 
                        PeriodicEntity ent(ids[other], o, local[other]);
                        m_periodicEdges[ids[i]].push_back(ent);
                    }
 
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        StdRegions::Orientation o =
                            orientMap[ids[i]] == orientMap[ids[other]]
                                ? StdRegions::eBackwards
                                : StdRegions::eForwards;
 
                        // Determine periodic vertices.
                        vector<int> perVerts1 = vertMap[ids[i]];
                        vector<int> perVerts2 = vertMap[ids[other]];
 
                        map<int, pair<int, bool>> tmpMap;
                        if (o == StdRegions::eForwards)
                        {
                            tmpMap[perVerts1[0]] = make_pair(
                                perVerts2[0], locVerts.count(perVerts2[0]) > 0);
                            tmpMap[perVerts1[1]] = make_pair(
                                perVerts2[1], locVerts.count(perVerts2[1]) > 0);
                        }
                        else
                        {
                            tmpMap[perVerts1[0]] = make_pair(
                                perVerts2[1], locVerts.count(perVerts2[1]) > 0);
                            tmpMap[perVerts1[1]] = make_pair(
                                perVerts2[0], locVerts.count(perVerts2[0]) > 0);
                        }
 
                        for (auto &mIt : tmpMap)
                        {
                            // See if this vertex has been recorded already.
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eNoOrientation,
                                                mIt.second.second);
                            auto perIt = periodicVerts.find(mIt.first);
 
                            if (perIt == periodicVerts.end())
                            {
                                periodicVerts[mIt.first].push_back(ent2);
                                perIt = periodicVerts.find(mIt.first);
                            }
                            else
                            {
                                bool doAdd = true;
                                for (j = 0; j < perIt->second.size(); ++j)
                                {
                                    if (perIt->second[j].id == mIt.second.first)
                                    {
                                        doAdd = false;
                                        break;
                                    }
                                }
 
                                if (doAdd)
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
                }
            }
 
            // Search for periodic vertices and edges which are not in
            // a periodic composite but lie in this process. First, loop
            // over all information we have from other processors.
            for (cnt = i = 0; i < totEdges; ++i)
            {
                int edgeId = edgeIds[i];
 
                ASSERTL0(allCompPairs.count(edgeId) > 0,
                         "Unable to find matching periodic edge.");
 
                int perEdgeId = allCompPairs[edgeId];
 
                for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
                {
                    int vId = traceIds[cnt];
 
                    auto perId = periodicVerts.find(vId);
 
                    if (perId == periodicVerts.end())
                    {
                        // This vertex is not included in the
                        // map. Figure out which vertex it is
                        // supposed to be periodic with. perEdgeId
                        // is the edge ID which is periodic with
                        // edgeId. The logic is much the same as
                        // the loop above.
                        int perVertexId =
                            orientMap[edgeId] == orientMap[perEdgeId]
                                ? vertMap[perEdgeId][(j + 1) % 2]
                                : vertMap[perEdgeId][j];
 
                        PeriodicEntity ent(perVertexId,
                                           StdRegions::eNoOrientation,
                                           locVerts.count(perVertexId) > 0);
 
                        periodicVerts[vId].push_back(ent);
                    }
                }
            }
 
            // Loop over all periodic vertices to complete connectivity
            // information.
            for (auto &perIt : periodicVerts)
            {
                // Loop over associated vertices.
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    auto perIt2 = periodicVerts.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicVerts.end(),
                             "Couldn't find periodic vertex.");
 
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        bool doAdd = true;
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                doAdd = false;
                                break;
                            }
                        }
 
                        if (doAdd)
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            // Do one final loop over periodic vertices to remove non-local
            // vertices from map.
            for (auto &perIt : periodicVerts)
            {
                if (locVerts.count(perIt.first) > 0)
                {
                    m_periodicVerts.insert(perIt);
                }
            }
        }
        break;
        case e3D:
        {
            SpatialDomains::CompositeOrdering compOrder =
                m_graph->GetCompositeOrdering();
            SpatialDomains::BndRegionOrdering bndRegOrder =
                m_graph->GetBndRegionOrdering();
            SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
 
            // perComps: Stores a unique collection of pairs of periodic
            // composites (i.e. if composites 1 and 2 are periodic then this map
            // will contain either the pair (1,2) or (2,1) but not both).
            //
            // The four maps allVerts, allCoord, allEdges and allOrient map a
            // periodic face to a vector containing the vertex ids of the face;
            // their coordinates; the edge ids of the face; and their
            // orientation within that face respectively.
            //
            // Finally the three sets locVerts, locEdges and locFaces store any
            // vertices, edges and faces that belong to a periodic composite and
            // lie on this process.
            map<int, RotPeriodicInfo> rotComp;
            map<int, int> perComps;
            map<int, vector<int>> allVerts;
            map<int, std::vector<SpatialDomains::PointGeom *>> allCoord;
            map<int, vector<int>> allEdges;
            map<int, vector<StdRegions::Orientation>> allOrient;
            set<int> locVerts;
            set<int> locEdges;
            set<int> locFaces;
 
            int region1ID, region2ID, i, j, k, cnt;
            SpatialDomains::BoundaryConditionShPtr locBCond;
 
            // Set up a set of all local verts and edges.
            for (i = 0; i < (*m_exp).size(); ++i)
            {
                for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetVid(j);
                    locVerts.insert(id);
                }
 
                for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetEid(j);
                    locEdges.insert(id);
                }
            }
 
            // Begin by populating the perComps map. We loop over all periodic
            // boundary conditions and determine the composite associated with
            // it, then fill out the all* maps.
            for (auto &it : bregions)
            {
 
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                // Check the region only contains a single composite.
                ASSERTL0(it.second->size() == 1,
                         "Boundary region " + std::to_string(region1ID) +
                             " should only contain 1 composite.");
 
                // From this identify composites by looking at the original
                // boundary region ordering. Note that in serial the mesh
                // partitioner is not run, so this map will be empty and
                // therefore needs to be populated by using the corresponding
                // boundary region.
                int cId1, cId2;
                if (vComm->GetSize() == 1)
                {
                    cId1 = it.second->begin()->first;
                    cId2 = bregions.find(region2ID)->second->begin()->first;
                }
                else
                {
                    cId1 = bndRegOrder.find(region1ID)->second[0];
                    cId2 = bndRegOrder.find(region2ID)->second[0];
                }
 
                // check to see if boundary is rotationally aligned
                if (boost::icontains(locBCond->GetUserDefined(), "Rotated"))
                {
                    vector<string> tmpstr;
 
                    boost::split(tmpstr, locBCond->GetUserDefined(),
                                 boost::is_any_of(":"));
 
                    if (boost::iequals(tmpstr[0], "Rotated"))
                    {
                        ASSERTL1(tmpstr.size() > 2,
                                 "Expected Rotated user defined string to "
                                 "contain direction and rotation angle "
                                 "and optionally a tolerance, "
                                 "i.e. Rotated:dir:PI/2:1e-6");
 
                        ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                                     (tmpstr[1] == "z"),
                                 "Rotated Dir is "
                                 "not specified as x,y or z");
 
                        RotPeriodicInfo RotInfo;
                        RotInfo.m_dir = (tmpstr[1] == "x")   ? 0
                                        : (tmpstr[1] == "y") ? 1
                                                             : 2;
 
                        LibUtilities::Interpreter strEval;
                        int ExprId      = strEval.DefineFunction("", tmpstr[2]);
                        RotInfo.m_angle = strEval.Evaluate(ExprId);
 
                        if (tmpstr.size() == 4)
                        {
                            try
                            {
                                RotInfo.m_tol = std::stod(tmpstr[3]);
                            }
                            catch (...)
                            {
                                NEKERROR(ErrorUtil::efatal,
                                         "failed to cast tolerance input "
                                         "to a double value in Rotated"
                                         "boundary information");
                            }
                        }
                        else
                        {
                            RotInfo.m_tol = 1e-8;
                        }
                        rotComp[cId1] = RotInfo;
                    }
                }
 
                SpatialDomains::CompositeSharedPtr c =
                    it.second->begin()->second;
 
                vector<unsigned int> tmpOrder;
 
                // store the rotation info of this
 
                // From the composite, we now construct the allVerts, allEdges
                // and allCoord map so that they can be transferred across
                // processors. We also populate the locFaces set to store a
                // record of all faces local to this process.
                for (i = 0; i < c->m_geomVec.size(); ++i)
                {
                    SpatialDomains::Geometry2D *faceGeom =
                        dynamic_cast<SpatialDomains::Geometry2D *>(
                            c->m_geomVec[i]);
                    ASSERTL1(faceGeom, "Unable to cast to Geometry2D*");
 
                    // Get geometry ID of this face and store in locFaces.
                    int faceId = c->m_geomVec[i]->GetGlobalID();
                    locFaces.insert(faceId);
 
                    // In serial, mesh partitioning will not have occurred so
                    // need to fill composite ordering map manually.
                    if (vComm->GetSize() == 1)
                    {
                        tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
                    }
 
                    // Loop over vertices and edges of the face to populate
                    // allVerts, allEdges and allCoord maps.
                    vector<int> vertList, edgeList;
                    std::vector<SpatialDomains::PointGeom *> coordVec;
                    vector<StdRegions::Orientation> orientVec;
                    for (j = 0; j < faceGeom->GetNumVerts(); ++j)
                    {
                        vertList.push_back(faceGeom->GetVid(j));
                        edgeList.push_back(faceGeom->GetEid(j));
                        coordVec.push_back(faceGeom->GetVertex(j));
                        orientVec.push_back(faceGeom->GetEorient(j));
                    }
 
                    allVerts[faceId]  = vertList;
                    allEdges[faceId]  = edgeList;
                    allCoord[faceId]  = coordVec;
                    allOrient[faceId] = orientVec;
                }
 
                // In serial, record the composite ordering in compOrder for
                // later in the routine.
                if (vComm->GetSize() == 1)
                {
                    compOrder[it.second->begin()->first] = tmpOrder;
                }
 
                // See if we already have either region1 or region2 stored in
                // perComps map already and do a sanity check to ensure regions
                // are mutually periodic.
                if (perComps.count(cId1) == 0)
                {
                    if (perComps.count(cId2) == 0)
                    {
                        perComps[cId1] = cId2;
                    }
                    else
                    {
                        std::stringstream ss;
                        ss << "Boundary region " << cId2 << " should be "
                           << "periodic with " << perComps[cId2] << " but "
                           << "found " << cId1 << " instead!";
                        ASSERTL0(perComps[cId2] == cId1, ss.str());
                    }
                }
                else
                {
                    std::stringstream ss;
                    ss << "Boundary region " << cId1 << " should be "
                       << "periodic with " << perComps[cId1] << " but "
                       << "found " << cId2 << " instead!";
                    ASSERTL0(perComps[cId1] == cId1, ss.str());
                }
            }
 
            // In case of periodic partition being split over many composites
            // may not have all periodic matches so check this
            int idmax = -1;
            for (auto &cIt : perComps)
            {
                idmax = max(idmax, cIt.first);
                idmax = max(idmax, cIt.second);
            }
            vComm->AllReduce(idmax, LibUtilities::ReduceMax);
            idmax++;
            Array<OneD, int> perid(idmax, -1);
            for (auto &cIt : perComps)
            {
                perid[cIt.first] = cIt.second;
            }
            vComm->AllReduce(perid, LibUtilities::ReduceMax);
            // update all partitions
            for (int i = 0; i < idmax; ++i)
            {
                if (perid[i] > -1)
                {
                    // skip if equivlaent relationship has
                    // already been speficied in list
                    if (perComps.count(perid[i]))
                    {
                        continue;
                    }
                    perComps[i] = perid[i];
                }
            }
 
            // The next routines process local face lists to
            // exchange vertices,
            // edges and faces.
            int n = vComm->GetSize();
            int p = vComm->GetRank();
            int totFaces;
            Array<OneD, int> facecounts(n, 0);
            Array<OneD, int> vertcounts(n, 0);
            Array<OneD, int> faceoffset(n, 0);
            Array<OneD, int> vertoffset(n, 0);
 
            Array<OneD, int> rotcounts(n, 0);
            Array<OneD, int> rotoffset(n, 0);
 
            rotcounts[p] = rotComp.size();
            vComm->AllReduce(rotcounts, LibUtilities::ReduceSum);
            int totrot = Vmath::Vsum(n, rotcounts, 1);
 
            if (totrot)
            {
                for (i = 1; i < n; ++i)
                {
                    rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
                }
 
                Array<OneD, int> compid(totrot, 0);
                Array<OneD, int> rotdir(totrot, 0);
                Array<OneD, NekDouble> rotangle(totrot, 0.0);
                Array<OneD, NekDouble> rottol(totrot, 0.0);
 
                // fill in rotational informaiton
                auto rIt = rotComp.begin();
 
                for (i = 0; rIt != rotComp.end(); ++rIt)
                {
                    compid[rotoffset[p] + i]   = rIt->first;
                    rotdir[rotoffset[p] + i]   = rIt->second.m_dir;
                    rotangle[rotoffset[p] + i] = rIt->second.m_angle;
                    rottol[rotoffset[p] + i++] = rIt->second.m_tol;
                }
 
                vComm->AllReduce(compid, LibUtilities::ReduceSum);
                vComm->AllReduce(rotdir, LibUtilities::ReduceSum);
                vComm->AllReduce(rotangle, LibUtilities::ReduceSum);
                vComm->AllReduce(rottol, LibUtilities::ReduceSum);
 
                // Fill in full rotational composite list
                for (i = 0; i < totrot; ++i)
                {
                    RotPeriodicInfo rinfo(rotdir[i], rotangle[i], rottol[i]);
 
                    rotComp[compid[i]] = rinfo;
                }
            }
 
            // First exchange the number of faces on each process.
            facecounts[p] = locFaces.size();
            vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
 
            // Set up an offset map to allow us to distribute face IDs to all
            // processors.
            faceoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
            }
 
            // Calculate total number of faces.
            totFaces = Vmath::Vsum(n, facecounts, 1);
 
            // faceIds holds face IDs for each periodic face. faceVerts holds
            // the number of vertices in this face.
            Array<OneD, int> faceIds(totFaces, 0);
            Array<OneD, int> faceVerts(totFaces, 0);
 
            // Process p writes IDs of its faces into position faceoffset[p] of
            // faceIds which allows us to perform an AllReduce to distribute
            // information amongst processors.
            auto sIt = locFaces.begin();
            for (i = 0; sIt != locFaces.end(); ++sIt)
            {
                faceIds[faceoffset[p] + i]     = *sIt;
                faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
            }
 
            vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
            vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
 
            // procVerts holds number of vertices (and also edges since each
            // face is 2D) on each process.
            Array<OneD, int> procVerts(n, 0);
            int nTotVerts;
 
            // Note if there are no periodic faces at all calling Vsum will
            // cause a segfault.
            if (totFaces > 0)
            {
                // Calculate number of vertices on each processor.
                nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
            }
            else
            {
                nTotVerts = 0;
            }
 
            for (i = 0; i < n; ++i)
            {
                if (facecounts[i] > 0)
                {
                    procVerts[i] = Vmath::Vsum(facecounts[i],
                                               faceVerts + faceoffset[i], 1);
                }
                else
                {
                    procVerts[i] = 0;
                }
            }
 
            // vertoffset is defined in the same manner as edgeoffset
            // beforehand.
            vertoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
            }
 
            // At this point we exchange all vertex IDs, edge IDs and vertex
            // coordinates for each face. The coordinates are necessary because
            // we need to calculate relative face orientations between periodic
            // faces to determined edge and vertex connectivity.
            Array<OneD, int> vertIds(nTotVerts, 0);
            Array<OneD, int> edgeIds(nTotVerts, 0);
            Array<OneD, int> edgeOrt(nTotVerts, 0);
            Array<OneD, NekDouble> vertX(nTotVerts, 0.0);
            Array<OneD, NekDouble> vertY(nTotVerts, 0.0);
            Array<OneD, NekDouble> vertZ(nTotVerts, 0.0);
 
            for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
            {
                for (j = 0; j < allVerts[*sIt].size(); ++j)
                {
                    int vertId                     = allVerts[*sIt][j];
                    vertIds[vertoffset[p] + cnt]   = vertId;
                    vertX[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(0);
                    vertY[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(1);
                    vertZ[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(2);
                    edgeIds[vertoffset[p] + cnt]   = allEdges[*sIt][j];
                    edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
                }
            }
 
            vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
            vComm->AllReduce(vertX, LibUtilities::ReduceSum);
            vComm->AllReduce(vertY, LibUtilities::ReduceSum);
            vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
 
            // Finally now we have all of this information, we construct maps
            // which make accessing the information easier. These are
            // conceptually the same as all* maps at the beginning of the
            // routine, but now hold information for all periodic vertices.
            map<int, vector<int>> vertMap;
            map<int, vector<int>> edgeMap;
            map<int, std::vector<SpatialDomains::PointGeom *>> coordMap;
 
            // These final two maps are required for determining the relative
            // orientation of periodic edges. vCoMap associates vertex IDs with
            // their coordinates, and eIdMap maps an edge ID to the two vertices
            // which construct it.
            map<int, SpatialDomains::PointGeomUniquePtr> vCoMap;
            map<int, pair<int, int>> eIdMap;
 
            for (cnt = i = 0; i < totFaces; ++i)
            {
                vector<int> edges(faceVerts[i]);
                vector<int> verts(faceVerts[i]);
                std::vector<SpatialDomains::PointGeom *> coord(faceVerts[i]);
 
                // Keep track of cnt to enable correct edge vertices to be
                // inserted into eIdMap.
                int tmp = cnt;
                for (j = 0; j < faceVerts[i]; ++j, ++cnt)
                {
                    edges[j] = edgeIds[cnt];
                    verts[j] = vertIds[cnt];
 
                    auto vIt = vCoMap.find(vertIds[cnt]);
 
                    if (vIt == vCoMap.end())
                    {
                        auto pt = ObjPoolManager<SpatialDomains::PointGeom>::
                            AllocateUniquePtr(3, verts[j], vertX[cnt],
                                              vertY[cnt], vertZ[cnt]);
                        coord[j]             = pt.get();
                        vCoMap[vertIds[cnt]] = std::move(pt);
                    }
                    else
                    {
                        coord[j] = vIt->second.get();
                    }
 
                    // Try to insert edge into the eIdMap to avoid re-inserting.
                    auto testIns = eIdMap.insert(make_pair(
                        edgeIds[cnt],
                        make_pair(vertIds[tmp + j],
                                  vertIds[tmp + ((j + 1) % faceVerts[i])])));
 
                    if (testIns.second == false)
                    {
                        continue;
                    }
 
                    // If the edge is reversed with respect to the face, then
                    // swap the edges so that we have the original ordering of
                    // the edge in the 3D element. This is necessary to properly
                    // determine edge orientation. Note that the logic relies on
                    // the fact that all edge forward directions are CCW
                    // orientated: we use a tensor product ordering for 2D
                    // elements so need to reverse this for edge IDs 2 and 3.
                    StdRegions::Orientation edgeOrient =
                        static_cast<StdRegions::Orientation>(edgeOrt[cnt]);
                    if (j > 1)
                    {
                        edgeOrient = edgeOrient == StdRegions::eForwards
                                         ? StdRegions::eBackwards
                                         : StdRegions::eForwards;
                    }
 
                    if (edgeOrient == StdRegions::eBackwards)
                    {
                        swap(testIns.first->second.first,
                             testIns.first->second.second);
                    }
                }
 
                vertMap[faceIds[i]]  = verts;
                edgeMap[faceIds[i]]  = edges;
                coordMap[faceIds[i]] = coord;
            }
 
            // Go through list of composites and figure out which edges are
            // parallel from original ordering in session file. This includes
            // composites which are not necessarily on this process.
 
            // Store temporary map of periodic vertices which will hold all
            // periodic vertices on the entire mesh so that doubly periodic
            // vertices/edges can be counted properly across partitions. Local
            // vertices/edges are copied into m_periodicVerts and
            // m_periodicEdges at the end of the function.
            PeriodicMap periodicVerts, periodicEdges;
 
            // Construct two maps which determine how vertices and edges of
            // faces connect given a specific face orientation. The key of the
            // map is the number of vertices in the face, used to determine
            // difference between tris and quads.
            map<int, map<StdRegions::Orientation, vector<int>>> vmap;
            map<int, map<StdRegions::Orientation, vector<int>>> emap;
 
            map<StdRegions::Orientation, vector<int>> quadVertMap;
            quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
            quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {3, 2, 1, 0};
            quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 3, 2};
            quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
            quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {0, 3, 2, 1};
            quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
            quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
            quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {2, 1, 0, 3};
 
            map<StdRegions::Orientation, vector<int>> quadEdgeMap;
            quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
            quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {2, 1, 0, 3};
            quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 3, 2, 1};
            quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
            quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {3, 2, 1, 0};
            quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
            quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
            quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {1, 0, 3, 2};
 
            map<StdRegions::Orientation, vector<int>> triVertMap;
            triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
            triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 2};
 
            map<StdRegions::Orientation, vector<int>> triEdgeMap;
            triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
            triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 2, 1};
 
            vmap[3] = triVertMap;
            vmap[4] = quadVertMap;
            emap[3] = triEdgeMap;
            emap[4] = quadEdgeMap;
 
            map<int, int> allCompPairs;
 
            // Collect composite id's of each periodic face for use if rotation
            // is required
            map<int, int> fIdToCompId;
 
            // Finally we have enough information to populate the periodic
            // vertex, edge and face maps. Begin by looping over all pairs of
            // periodic composites to determine pairs of periodic faces.
            for (auto &cIt : perComps)
            {
                SpatialDomains::CompositeSharedPtr c[2];
                const int id1    = cIt.first;
                const int id2    = cIt.second;
                std::string id1s = std::to_string(id1);
                std::string id2s = std::to_string(id2);
 
                if (compMap.count(id1) > 0)
                {
                    c[0] = compMap[id1];
                }
 
                if (compMap.count(id2) > 0)
                {
                    c[1] = compMap[id2];
                }
 
                // Loop over composite ordering to construct list of all
                // periodic faces, regardless of whether they are on this
                // process.
                map<int, int> compPairs;
 
                ASSERTL0(compOrder.count(id1) > 0,
                         "Unable to find composite " + id1s + " in order map.");
                ASSERTL0(compOrder.count(id2) > 0,
                         "Unable to find composite " + id2s + " in order map.");
                WARNINGL0(
                    compOrder[id1].size() == compOrder[id2].size(),
                    "Periodic composites " + id1s + " and " + id2s +
                        " should have the same number of elements. Have " +
                        std::to_string(compOrder[id1].size()) + " vs " +
                        std::to_string(compOrder[id2].size()));
                ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
                                                        id1s + " and " + id2s +
                                                        " are empty!");
 
                // Look up composite ordering to determine pairs.
                for (i = 0; i < compOrder[id1].size(); ++i)
                {
                    int eId1 = compOrder[id1][i];
                    int eId2 = compOrder[id2][i];
 
                    ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
 
                    // Sanity check that the faces are mutually periodic.
                    if (compPairs.count(eId2) != 0)
                    {
                        ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
                    }
                    compPairs[eId1] = eId2;
 
                    // store  a map of face ids to composite ids
                    fIdToCompId[eId1] = id1;
                    fIdToCompId[eId2] = id2;
                }
 
                // Now that we have all pairs of periodic faces, loop over the
                // ones local on this process and populate face/edge/vertex
                // maps.
                for (auto &pIt : compPairs)
                {
                    int ids[2]    = {pIt.first, pIt.second};
                    bool local[2] = {locFaces.count(pIt.first) > 0,
                                     locFaces.count(pIt.second) > 0};
 
                    ASSERTL0(coordMap.count(ids[0]) > 0 &&
                                 coordMap.count(ids[1]) > 0,
                             "Unable to find face in coordinate map");
 
                    allCompPairs[pIt.first]  = pIt.second;
                    allCompPairs[pIt.second] = pIt.first;
 
                    // Loop up coordinates of the faces, check they have the
                    // same number of vertices.
                    std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                        coordMap[ids[0]], coordMap[ids[1]]};
 
                    ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
                             "Two periodic faces have different number "
                             "of vertices!");
 
                    // o will store relative orientation of faces. Note that in
                    // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
                    // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
                    // different going from face1->face2 instead of face2->face1
                    // (check this).
                    StdRegions::Orientation o;
                    bool rotbnd     = false;
                    int dir         = 0;
                    NekDouble angle = 0.0;
                    NekDouble sign  = 0.0;
                    NekDouble tol   = 1e-8;
 
                    // check to see if perioid boundary is rotated
                    if (rotComp.count(fIdToCompId[pIt.first]))
                    {
                        rotbnd = true;
                        dir    = rotComp[fIdToCompId[pIt.first]].m_dir;
                        angle  = rotComp[fIdToCompId[pIt.first]].m_angle;
                        tol    = rotComp[fIdToCompId[pIt.first]].m_tol;
                    }
 
                    // Record periodic faces.
                    for (i = 0; i < 2; ++i)
                    {
                        if (!local[i])
                        {
                            continue;
                        }
 
                        // Reference to the other face.
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        // Calculate relative face orientation.
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        // Record face ID, orientation and whether other face is
                        // local.
                        PeriodicEntity ent(ids[other], o, local[other]);
                        m_periodicFaces[ids[i]].push_back(ent);
                    }
 
                    int nFaceVerts = vertMap[ids[0]].size();
 
                    // Determine periodic vertices.
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        // Calculate relative face orientation.
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        if (nFaceVerts == 3)
                        {
                            ASSERTL0(
                                o == StdRegions::eDir1FwdDir1_Dir2FwdDir2 ||
                                    o == StdRegions::eDir1BwdDir1_Dir2FwdDir2,
                                "Unsupported face orientation for face " +
                                    std::to_string(ids[i]));
                        }
 
                        // Look up vertices for this face.
                        vector<int> per1 = vertMap[ids[i]];
                        vector<int> per2 = vertMap[ids[other]];
 
                        // tmpMap will hold the pairs of vertices which are
                        // periodic.
                        map<int, pair<int, bool>> tmpMap;
 
                        // Use vmap to determine which vertices connect given
                        // the orientation o.
                        for (j = 0; j < nFaceVerts; ++j)
                        {
                            int v = vmap[nFaceVerts][o][j];
                            tmpMap[per1[j]] =
                                make_pair(per2[v], locVerts.count(per2[v]) > 0);
                        }
 
                        // Now loop over tmpMap to associate periodic vertices.
                        for (auto &mIt : tmpMap)
                        {
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eNoOrientation,
                                                mIt.second.second);
 
                            // See if this vertex has been recorded already.
                            auto perIt = periodicVerts.find(mIt.first);
 
                            if (perIt == periodicVerts.end())
                            {
                                // Vertex is new - just record this entity as
                                // usual.
                                periodicVerts[mIt.first].push_back(ent2);
                                perIt = periodicVerts.find(mIt.first);
                            }
                            else
                            {
                                // Vertex is known - loop over the vertices
                                // inside the record and potentially add vertex
                                // mIt.second to the list.
                                for (k = 0; k < perIt->second.size(); ++k)
                                {
                                    if (perIt->second[k].id == mIt.second.first)
                                    {
                                        break;
                                    }
                                }
 
                                if (k == perIt->second.size())
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
 
                    // Determine periodic edges. Logic is the same as above,
                    // and perhaps should be condensed to avoid replication.
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        vector<int> per1 = edgeMap[ids[i]];
                        vector<int> per2 = edgeMap[ids[other]];
 
                        map<int, pair<int, bool>> tmpMap;
 
                        for (j = 0; j < nFaceVerts; ++j)
                        {
                            int e = emap[nFaceVerts][o][j];
                            tmpMap[per1[j]] =
                                make_pair(per2[e], locEdges.count(per2[e]) > 0);
                        }
 
                        for (auto &mIt : tmpMap)
                        {
                            // Note we assume orientation of edges is forwards -
                            // this may be reversed later.
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eForwards,
                                                mIt.second.second);
                            auto perIt = periodicEdges.find(mIt.first);
 
                            if (perIt == periodicEdges.end())
                            {
                                periodicEdges[mIt.first].push_back(ent2);
                                perIt = periodicEdges.find(mIt.first);
                            }
                            else
                            {
                                for (k = 0; k < perIt->second.size(); ++k)
                                {
                                    if (perIt->second[k].id == mIt.second.first)
                                    {
                                        break;
                                    }
                                }
 
                                if (k == perIt->second.size())
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
                }
            }
 
            // Make sure that the nubmer of face pairs and the
            // face Id to composite Id map match in size
            ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
                     "At this point the size of allCompPairs "
                     "should have been the same as fIdToCompId");
 
            // also will need an edge id to composite id at
            // end of routine
            map<int, int> eIdToCompId;
 
            // Search for periodic vertices and edges which are not
            // in a periodic composite but lie in this process. First,
            // loop over all information we have from other
            // processors.
            for (cnt = i = 0; i < totFaces; ++i)
            {
                bool rotbnd     = false;
                int dir         = 0;
                NekDouble angle = 0.0;
                NekDouble tol   = 1e-8;
 
                int faceId = faceIds[i];
 
                ASSERTL0(allCompPairs.count(faceId) > 0,
                         "Unable to find matching periodic face. faceId = " +
                             std::to_string(faceId));
 
                int perFaceId = allCompPairs[faceId];
 
                // check to see if periodic boundary is rotated
                ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
                         "Face " + std::to_string(faceId) +
                             " not found in fIdtoCompId map");
                if (rotComp.count(fIdToCompId[faceId]))
                {
                    rotbnd = true;
                    dir    = rotComp[fIdToCompId[faceId]].m_dir;
                    angle  = rotComp[fIdToCompId[faceId]].m_angle;
                    tol    = rotComp[fIdToCompId[faceId]].m_tol;
                }
 
                for (j = 0; j < faceVerts[i]; ++j, ++cnt)
                {
                    int vId = vertIds[cnt];
 
                    auto perId = periodicVerts.find(vId);
 
                    if (perId == periodicVerts.end())
                    {
 
                        // This vertex is not included in the
                        // map. Figure out which vertex it is supposed
                        // to be periodic with. perFaceId is the face
                        // ID which is periodic with faceId. The logic
                        // is much the same as the loop above.
                        std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                            coordMap[faceId], coordMap[perFaceId]};
 
                        int nFaceVerts = tmpVec[0].size();
                        StdRegions::Orientation o =
                            nFaceVerts == 3
                                ? SpatialDomains::TriGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1],
                                       tmpVec[0][2]},
                                      {tmpVec[1][0], tmpVec[1][1],
                                       tmpVec[1][2]},
                                      rotbnd, dir, angle, tol)
                                : SpatialDomains::QuadGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1], tmpVec[0][2],
                                       tmpVec[0][3]},
                                      {tmpVec[1][0], tmpVec[1][1], tmpVec[1][2],
                                       tmpVec[1][3]},
                                      rotbnd, dir, angle, tol);
 
                        // Use vmap to determine which vertex of the other face
                        // should be periodic with this one.
                        int perVertexId =
                            vertMap[perFaceId][vmap[nFaceVerts][o][j]];
 
                        PeriodicEntity ent(perVertexId,
                                           StdRegions::eNoOrientation,
                                           locVerts.count(perVertexId) > 0);
 
                        periodicVerts[vId].push_back(ent);
                    }
 
                    int eId = edgeIds[cnt];
 
                    perId = periodicEdges.find(eId);
 
                    // this map is required at very end to determine rotation of
                    // edges.
                    if (rotbnd)
                    {
                        eIdToCompId[eId] = fIdToCompId[faceId];
                    }
 
                    if (perId == periodicEdges.end())
                    {
                        // This edge is not included in the map. Figure
                        // out which edge it is supposed to be periodic
                        // with. perFaceId is the face ID which is
                        // periodic with faceId. The logic is much the
                        // same as the loop above.
                        std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                            coordMap[faceId], coordMap[perFaceId]};
 
                        int nFaceEdges = tmpVec[0].size();
                        StdRegions::Orientation o =
                            nFaceEdges == 3
                                ? SpatialDomains::TriGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1],
                                       tmpVec[0][2]},
                                      {tmpVec[1][0], tmpVec[1][1],
                                       tmpVec[1][2]},
                                      rotbnd, dir, angle, tol)
                                : SpatialDomains::QuadGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1], tmpVec[0][2],
                                       tmpVec[0][3]},
                                      {tmpVec[1][0], tmpVec[1][1], tmpVec[1][2],
                                       tmpVec[1][3]},
                                      rotbnd, dir, angle, tol);
 
                        // Use emap to determine which edge of the other
                        // face should be periodic with this one.
                        int perEdgeId =
                            edgeMap[perFaceId][emap[nFaceEdges][o][j]];
 
                        PeriodicEntity ent(perEdgeId, StdRegions::eForwards,
                                           locEdges.count(perEdgeId) > 0);
 
                        periodicEdges[eId].push_back(ent);
 
                        // this map is required at very end to
                        // determine rotation of edges.
                        if (rotbnd)
                        {
                            eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
                        }
                    }
                }
            }
 
            // Finally, we must loop over the periodicVerts and periodicEdges
            // map to complete connectivity information.
            for (auto &perIt : periodicVerts)
            {
                // For each vertex that is periodic with this one...
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    // Find the vertex in the periodicVerts map...
                    auto perIt2 = periodicVerts.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicVerts.end(),
                             "Couldn't find periodic vertex.");
 
                    // Now search through this vertex's list and make sure that
                    // we have a record of any vertices which aren't in the
                    // original list.
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                break;
                            }
                        }
 
                        if (k == perIt.second.size())
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            for (auto &perIt : periodicEdges)
            {
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    auto perIt2 = periodicEdges.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicEdges.end(),
                             "Couldn't find periodic edge.");
 
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                break;
                            }
                        }
 
                        if (k == perIt.second.size())
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            // Loop over periodic edges to determine relative edge orientations.
            for (auto &perIt : periodicEdges)
            {
                bool rotbnd     = false;
                int dir         = 0;
                NekDouble angle = 0.0;
                NekDouble tol   = 1e-8;
 
                // Find edge coordinates
                auto eIt                       = eIdMap.find(perIt.first);
                SpatialDomains::PointGeom v[2] = {*vCoMap[eIt->second.first],
                                                  *vCoMap[eIt->second.second]};
 
                // check to see if perioid boundary is rotated
                if (rotComp.count(eIdToCompId[perIt.first]))
                {
                    rotbnd = true;
                    dir    = rotComp[eIdToCompId[perIt.first]].m_dir;
                    angle  = rotComp[eIdToCompId[perIt.first]].m_angle;
                    tol    = rotComp[eIdToCompId[perIt.first]].m_tol;
                }
 
                // Loop over each edge, and construct a vector that takes us
                // from one vertex to another. Use this to figure out which
                // vertex maps to which.
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    eIt = eIdMap.find(perIt.second[i].id);
 
                    SpatialDomains::PointGeom w[2] = {
                        *vCoMap[eIt->second.first],
                        *vCoMap[eIt->second.second]};
 
                    int vMap[2] = {-1, -1};
                    if (rotbnd)
                    {
 
                        SpatialDomains::PointGeom r;
 
                        r.Rotate(v[0], dir, angle);
 
                        if (r.dist(w[0]) < tol)
                        {
                            vMap[0] = 0;
                        }
                        else
                        {
                            r.Rotate(v[1], dir, angle);
                            if (r.dist(w[0]) < tol)
                            {
                                vMap[0] = 1;
                            }
                            else
                            {
                                NEKERROR(ErrorUtil::efatal,
                                         "Unable to align rotationally "
                                         "periodic edge vertex");
                            }
                        }
                    }
                    else // translation test
                    {
                        NekDouble cx =
                            0.5 * (w[0](0) - v[0](0) + w[1](0) - v[1](0));
                        NekDouble cy =
                            0.5 * (w[0](1) - v[0](1) + w[1](1) - v[1](1));
                        NekDouble cz =
                            0.5 * (w[0](2) - v[0](2) + w[1](2) - v[1](2));
 
                        for (j = 0; j < 2; ++j)
                        {
                            NekDouble x = v[j](0);
                            NekDouble y = v[j](1);
                            NekDouble z = v[j](2);
                            for (k = 0; k < 2; ++k)
                            {
                                NekDouble x1 = w[k](0) - cx;
                                NekDouble y1 = w[k](1) - cy;
                                NekDouble z1 = w[k](2) - cz;
 
                                if (sqrt((x1 - x) * (x1 - x) +
                                         (y1 - y) * (y1 - y) +
                                         (z1 - z) * (z1 - z)) < 1e-8)
                                {
                                    vMap[k] = j;
                                    break;
                                }
                            }
                        }
 
                        // Sanity check the map.
                        ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
                                 "Unable to align periodic edge vertex.");
                        ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
                                     (vMap[1] == 0 || vMap[1] == 1) &&
                                     (vMap[0] != vMap[1]),
                                 "Unable to align periodic edge vertex.");
                    }
 
                    // If 0 -> 0 then edges are aligned already; otherwise
                    // reverse the orientation.
                    if (vMap[0] != 0)
                    {
                        perIt.second[i].orient = StdRegions::eBackwards;
                    }
                }
            }
 
            // Do one final loop over periodic vertices/edges to remove
            // non-local vertices/edges from map.
            for (auto &perIt : periodicVerts)
            {
                if (locVerts.count(perIt.first) > 0)
                {
                    m_periodicVerts.insert(perIt);
                }
            }
 
            for (auto &perIt : periodicEdges)
            {
                if (locEdges.count(perIt.first) > 0)
                {
                    m_periodicEdges.insert(perIt);
                }
            }
        }
        break;
        default:
            ASSERTL1(false, "Not setup for this expansion");
            break;
    }
}

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), 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::Geometry::GetEid(), Nektar::SpatialDomains::Geometry::GetEorient(), Nektar::SpatialDomains::QuadGeom::GetFaceOrientation(), Nektar::SpatialDomains::TriGeom::GetFaceOrientation(), Nektar::SpatialDomains::Geometry::GetNumVerts(), Nektar::SpatialDomains::Geometry::GetVertex(), Nektar::SpatialDomains::Geometry::GetVid(), 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, tinysimd::max(), NEKERROR, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Nektar::SpatialDomains::PointGeom::Rotate(), sign, tinysimd::sqrt(), Vmath::Vsum(), and WARNINGL0.

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

GenerateBoundaryConditionExpansion() [1/2]

void Nektar::MultiRegions::DisContField::GenerateBoundaryConditionExpansion ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  In, const SpatialDomains::BoundaryConditions &  bcs, const std::string  variable, const bool  DeclareCoeffPhysArrays = true, const Collections::ImplementationType  ImpType = Collections::eNoImpType  )

protected

Make copy of boundary conditions.

This function discretises the boundary conditions using an already setup expansion 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

In	An array of boundary conditions from an alread setup expansion and the spectral/hp element expansions.
bcs	Information about the enforced boundary conditions.
variable	The session variable associated with the boundary conditions to enforce.
DeclareCoeffPhysArrays	bool to identify if array space should be setup. Default is true.

Definition at line 918 of file DisContField.cpp.

{
    int cnt = 0;
    SpatialDomains::BoundaryConditionShPtr bc;
    MultiRegions::ExpListSharedPtr locExpList;
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    std::set<int> ProcessBnd;
    if (m_session->DefinesTag("CreateBndRegions"))
    {
        // evaluate bnd regions to be generated
        vector<unsigned int> bndRegions;
        ASSERTL0(ParseUtils::GenerateVector(
                     m_session->GetTag("CreateBndRegions"), bndRegions),
                 "Failed to interpret bnd values string");
 
        for (auto &bnd : bndRegions)
        {
            if (bregions.count(bnd) == 1)
            {
                ProcessBnd.insert(bnd);
            }
        }
    }
    else
    {
        for (auto &it : bregions)
        {
            ProcessBnd.insert(it.first);
        }
    }
 
    m_bndCondExpansions =
        Array<OneD, MultiRegions::ExpListSharedPtr>(ProcessBnd.size());
    m_bndConditions =
        Array<OneD, SpatialDomains::BoundaryConditionShPtr>(ProcessBnd.size());
 
    m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
 
    // count the number of non-periodic boundary points
    for (auto &it : bregions)
    {
        // check to see if reduced bnd regions have been set in FieldConvert.
        if (ProcessBnd.count(it.first) == 0)
        {
            continue;
        }
 
        bc = GetBoundaryCondition(bconditions, it.first, variable);
 
        // make a copy of existing Bc passed to generate function
        locExpList = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
            *(In[cnt]), DeclareCoeffPhysArrays);
 
        m_bndCondExpansions[cnt] = locExpList;
        m_bndConditions[cnt]     = bc;
 
        std::string type = m_bndConditions[cnt]->GetUserDefined();
 
        // Set up normals on non-Dirichlet boundary conditions. Second
        // two conditions ideally should be in local solver setup (when
        // made into factory)
        if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
            boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
        {
            SetUpPhysNormals();
        }
        cnt++;
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::ParseUtils::GenerateVector(), 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().

GenerateBoundaryConditionExpansion() [2/2]

void Nektar::MultiRegions::DisContField::GenerateBoundaryConditionExpansion ( const SpatialDomains::MeshGraphSharedPtr &  graph, const SpatialDomains::BoundaryConditions &  bcs, const std::string  variable, const bool  DeclareCoeffPhysArrays = true, const Collections::ImplementationType  ImpType = Collections::eNoImpType  )

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

graph	A mesh, containing information about the domain and the spectral/hp element expansions.
bcs	Information about the enforced boundary conditions.
variable	The session variable associated with the boundary conditions to enforce.
DeclareCoeffPhysArrays	bool to identify if array space should be setup. Default is true.

Definition at line 822 of file DisContField.cpp.

{
    int cnt = 0;
    SpatialDomains::BoundaryConditionShPtr bc;
    MultiRegions::ExpListSharedPtr locExpList;
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    std::set<int> ProcessBnd;
    if (m_session->DefinesTag("CreateBndRegions"))
    {
        // evaluate bnd regions to be generated
        vector<unsigned int> bndRegions;
        ASSERTL0(ParseUtils::GenerateVector(
                     m_session->GetTag("CreateBndRegions"), bndRegions),
                 "Failed to interpret bnd values string");
 
        for (auto &bnd : bndRegions)
        {
            if (bregions.count(bnd) ==
                1) // this boundary may not exist on processor
            {
                ProcessBnd.insert(bnd);
            }
        }
    }
    else
    {
        for (auto &it : bregions)
        {
            ProcessBnd.insert(it.first);
        }
    }
 
    m_bndCondExpansions =
        Array<OneD, MultiRegions::ExpListSharedPtr>(ProcessBnd.size());
    m_bndConditions =
        Array<OneD, SpatialDomains::BoundaryConditionShPtr>(ProcessBnd.size());
 
    m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
 
    // count the number of non-periodic boundary points
    for (auto &it : bregions)
    {
        // check to see if reduced bnd regions have been set in FieldConvert.
        if (ProcessBnd.count(it.first) == 0)
        {
            continue;
        }
 
        bc = GetBoundaryCondition(bconditions, it.first, variable);
 
        locExpList = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
            m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
            false, bc->GetComm(), ImpType);
 
        m_bndCondExpansions[cnt] = locExpList;
        m_bndConditions[cnt]     = bc;
 
        std::string type = m_bndConditions[cnt]->GetUserDefined();
 
        // Set up normals on non-Dirichlet boundary conditions. Second
        // two conditions ideally should be in local solver setup (when
        // made into factory)
        if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
            boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
        {
            SetUpPhysNormals();
        }
        cnt++;
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::ParseUtils::GenerateVector(), 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(), DisContField(), and 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 518 of file DisContField.cpp.

{
    SpatialDomains::BoundaryConditionsSharedPtr returnval;
 
    returnval =
        MemoryManager<SpatialDomains::BoundaryConditions>::AllocateSharedPtr();
 
    map<int, int> GeometryToRegionsMap;
 
    const SpatialDomains::BoundaryRegionCollection &bregions =
        Allbcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        Allbcs.GetBoundaryConditions();
 
    // Set up a map of all boundary regions
    for (auto &it : bregions)
    {
        for (auto &bregionIt : *it.second)
        {
            // can assume that all regions only contain one point in 1D
            // Really do not need loop above
            int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
            GeometryToRegionsMap[id] = it.first;
        }
    }
 
    map<int, SpatialDomains::Geometry *> EndOfDomain;
 
    // Now find out which points in domain have only one vertex
    for (auto &domIt : domain)
    {
        SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
        for (int i = 0; i < geomvector->m_geomVec.size(); ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                int vid = geomvector->m_geomVec[i]->GetVid(j);
                if (EndOfDomain.count(vid) == 0)
                {
                    EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
                }
                else
                {
                    EndOfDomain.erase(vid);
                }
            }
        }
    }
    ASSERTL1(EndOfDomain.size() == 2, "Did not find two ends of domain");
 
    for (auto &regIt : EndOfDomain)
    {
        if (GeometryToRegionsMap.count(regIt.first) != 0)
        {
            // Set up boundary condition up
            auto iter = GeometryToRegionsMap.find(regIt.first);
            ASSERTL1(iter != GeometryToRegionsMap.end(),
                     "Failied to find GeometryToRegionMap");
 
            int regionId      = iter->second;
            auto bregionsIter = bregions.find(regionId);
            ASSERTL1(bregionsIter != bregions.end(),
                     "Failed to find boundary region");
 
            SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
            returnval->AddBoundaryRegions(regionId, breg);
 
            auto bconditionsIter = bconditions.find(regionId);
            ASSERTL1(bconditionsIter != bconditions.end(),
                     "Failed to find boundary collection");
 
            SpatialDomains::BoundaryConditionMapShPtr bcond =
                bconditionsIter->second;
            returnval->AddBoundaryConditions(regionId, bcond);
        }
    }
 
    return returnval;
}

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

field

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

{
    ASSERTL1((*m_exp)[0]->GetBasis(0)->GetBasisType() !=
                 LibUtilities::eGauss_Lagrange,
             "Routine not set up for Gauss Lagrange points distribution");
 
    // blocked routine
    Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->GetNLocTracePts());
 
    m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
    m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
 
    Array<OneD, NekDouble> invals =
        tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
    m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
 
    DisContField::v_PeriodicBwdCopy(Fwd, Bwd);
 
    FillBwdWithBoundCond(Fwd, Bwd, BndCond, BndCondExp, false);
 
    // Do parallel exchange for forwards/backwards spaces.
    m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
 
    // Do exchange of interface traces (local and parallel)
    // We may have to split this out into separate local and
    // parallel for IP method???
    m_interfaceMap->ExchangeTrace(Fwd, Bwd);
}

References ASSERTL1, Nektar::LibUtilities::eGauss_Lagrange, 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 2886 of file DisContField.cpp.

{
    ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
             "Routine currently only tested for HybridDGHelmholtz");
 
    ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
             "Full matrix global systems are not supported for HDG "
             "expansions");
 
    ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
             "The local to global map is not set up for the requested "
             "solution type");
 
    GlobalLinSysSharedPtr glo_matrix;
    auto matrixIter = m_globalBndMat->find(mkey);
 
    if (matrixIter == m_globalBndMat->end())
    {
        glo_matrix              = GenGlobalBndLinSys(mkey, m_traceMap);
        (*m_globalBndMat)[mkey] = glo_matrix;
    }
    else
    {
        glo_matrix = matrixIter->second;
    }
 
    return glo_matrix;
}

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

GetLocElmtTrace()

ExpListSharedPtr & Nektar::MultiRegions::DisContField::GetLocElmtTrace ( )

inline

Definition at line 142 of file DisContField.h.

    {
        if (!m_locElmtTrace)
        {
            m_locElmtTrace =
                MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
                    m_session, *m_exp, GetGraph(), true, "LocElmtTrace");
        }
 
        return m_locElmtTrace;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::ExpList::GetGraph(), Nektar::MultiRegions::ExpList::m_exp, m_locElmtTrace, and Nektar::MultiRegions::ExpList::m_session.

Referenced by Nektar::SolverUtils::DiffusionIP::AddSymmFluxIntegralToCoeff().

GetLocTraceToTraceMap()

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

inline

Definition at line 123 of file DisContField.h.

    {
        LocTraceToTraceMap = m_locTraceToTraceMap;
    }

References m_locTraceToTraceMap.

GetNegatedFluxNormal()

vector< bool > & Nektar::MultiRegions::DisContField::GetNegatedFluxNormal

(

void

)

Definition at line 2950 of file DisContField.cpp.

{
    if (m_negatedFluxNormal.size() == 0)
    {
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
 
        m_negatedFluxNormal.resize(2 * GetExpSize());
 
        for (int i = 0; i < GetExpSize(); ++i)
        {
 
            for (int v = 0; v < 2; ++v)
            {
                LocalRegions::Expansion0DSharedPtr vertExp =
                    elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
 
                if (vertExp->GetLeftAdjacentElementExp()
                        ->GetGeom()
                        ->GetGlobalID() !=
                    (*m_exp)[i]->GetGeom()->GetGlobalID())
                {
                    m_negatedFluxNormal[2 * i + v] = true;
                }
                else
                {
                    m_negatedFluxNormal[2 * i + v] = false;
                }
            }
        }
    }
 
    return m_negatedFluxNormal;
}

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

Referenced by v_AddTraceIntegral().

GetTraceElmtId()

unsigned Nektar::MultiRegions::DisContField::GetTraceElmtId ( const unsigned  elmtid, const unsigned  traceid  )

inline

Definition at line 136 of file DisContField.h.

    {
        return m_traceMap->GetElmtToTrace()[elmtid][traceid]->GetElmtId();
    }

References m_traceMap.

IsLeftAdjacentTrace()

bool Nektar::MultiRegions::DisContField::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.

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

{
    LocalRegions::ExpansionSharedPtr traceEl =
        m_traceMap->GetElmtToTrace()[n][e];
 
    PeriodicMap periodicTraces = (m_expType == e1D)   ? m_periodicVerts
                                 : (m_expType == e2D) ? m_periodicEdges
                                                      : m_periodicFaces;
 
    bool fwd = true;
    if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
        traceEl->GetRightAdjacentElementTrace() == -1)
    {
        // Boundary edge (1 connected element). Do nothing in
        // serial.
        auto it = m_boundaryTraces.find(traceEl->GetElmtId());
 
        // If the edge does not have a boundary condition set on
        // it, then assume it is a partition edge or periodic.
        if (it == m_boundaryTraces.end())
        {
            fwd = true;
        }
    }
    else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
             traceEl->GetRightAdjacentElementTrace() != -1)
    {
        // Non-boundary edge (2 connected elements).
        fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
    }
    else
    {
        ASSERTL2(false, "Unconnected trace element!");
    }
 
    return fwd;
}

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

{
    int i, e, ncoeff_edge;
    Array<OneD, const NekDouble> tmp_coeffs;
    Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
 
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    StdRegions::Orientation edgedir;
 
    int cnt;
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
    m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
 
    edge_lambda = loc_lambda;
 
    // Calculate Q using standard DG formulation.
    for (i = cnt = 0; i < GetExpSize(); ++i)
    {
        // Probably a better way of setting up lambda than this.
        // Note cannot use PutCoeffsInToElmts since lambda space
        // is mapped during the solve.
        int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
        Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
 
        for (e = 0; e < nEdges; ++e)
        {
            edgedir       = (*m_exp)[i]->GetTraceOrient(e);
            ncoeff_edge   = elmtToTrace[i][e]->GetNcoeffs();
            edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
            Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
            elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
                                                      edgeCoeffs[e]);
            edge_lambda = edge_lambda + ncoeff_edge;
        }
 
        (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
        cnt += (*m_exp)[i]->GetNcoeffs();
    }
 
    Array<OneD, NekDouble> phys(m_npoints);
    BwdTrans(out_d, phys);
    Vmath::Vsub(m_npoints, phys, 1, soln, 1, phys, 1);
    return L2(phys);
}

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

Rotate()

void Nektar::MultiRegions::DisContField::Rotate ( Array< OneD, Array< OneD, NekDouble > > &  Bwd, const int  dir, const NekDouble  angle, const int  offset, const int  npts  )

protected

Rotate the slice [offset, offset + npts) of the 3D vector field in Bwd around the axis 'dir' by 'angle'.

Definition at line 4877 of file DisContField.cpp.

{
    // Precompute trigonometric values for improved precision
    const NekDouble cosA = cos(angle);
    const NekDouble sinA = sin(angle);
 
    switch (dir)
    {
        // ----------------------------------------------------------
        // Rotate around the x-axis by 'angle':
        //  y' =  y*cos(angle) - z*sin(angle)
        //  z' =  y*sin(angle) + z*cos(angle)
        //  x' =  x (unchanged)
        // ----------------------------------------------------------
        case 0:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpY = Bwd[2][i];
                NekDouble tmpZ = Bwd[3][i];
 
                NekDouble yrot = cosA * tmpY - sinA * tmpZ;
                NekDouble zrot = sinA * tmpY + cosA * tmpZ;
 
                Bwd[2][i] = yrot;
                Bwd[3][i] = zrot;
            }
        }
        break;
 
        // ----------------------------------------------------------
        // Rotate around the y-axis by 'angle':
        //  z' =  z*cos(angle) - x*sin(angle)
        //  x' =  z*sin(angle) + x*cos(angle)
        //  y' =  y (unchanged)
        // ----------------------------------------------------------
        case 1:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpX = Bwd[1][i];
                NekDouble tmpZ = Bwd[3][i];
 
                NekDouble zrot = cosA * tmpZ - sinA * tmpX;
                NekDouble xrot = sinA * tmpZ + cosA * tmpX;
 
                Bwd[1][i] = xrot;
                Bwd[3][i] = zrot;
            }
        }
        break;
 
        // ----------------------------------------------------------
        // Rotate around the z-axis by 'angle':
        //  x' =  x*cos(angle) - y*sin(angle)
        //  y' =  x*sin(angle) + y*cos(angle)
        //  z' =  z (unchanged)
        // ----------------------------------------------------------
        case 2:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpX = Bwd[1][i];
                NekDouble tmpY = Bwd[2][i];
 
                NekDouble xrot = cosA * tmpX - sinA * tmpY;
                NekDouble yrot = sinA * tmpX + cosA * tmpY;
 
                Bwd[1][i] = xrot;
                Bwd[2][i] = yrot;
            }
        }
        break;
 
        default:
            NEKERROR(ErrorUtil::efatal,
                     "Rotate() axis must be 0 (x), 1 (y), or 2 (z).");
            break;
    }
}

References Nektar::ErrorUtil::efatal, and NEKERROR.

Referenced by v_PeriodicBwdRot().

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

In	Field to compare with.

Returns

True if boundary conditions match.

Definition at line 2922 of file DisContField.cpp.

{
    int i;
    bool returnval = true;
 
    for (i = 0; i < m_bndConditions.size(); ++i)
    {
        // check to see if boundary condition type is the same
        // and there are the same number of boundary
        // conditions in the boundary definition.
        if ((m_bndConditions[i]->GetBoundaryConditionType() !=
             In.m_bndConditions[i]->GetBoundaryConditionType()) ||
            (m_bndCondExpansions[i]->GetExpSize() !=
             In.m_bndCondExpansions[i]->GetExpSize()))
        {
            returnval = false;
            break;
        }
    }
 
    // Compare with all other processes. Return true only if all
    // processes report having the same boundary conditions.
    int vSame = (returnval ? 1 : 0);
    m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
 
    return (vSame == 1);
}

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

{
    // Check for multiple calls
    if (m_trace != NullExpListSharedPtr)
    {
        return;
    }
 
    // Set up matrix map
    m_globalBndMat = MemoryManager<GlobalLinSysMap>::AllocateSharedPtr();
 
    // Set up trace space
    m_trace = MemoryManager<ExpList>::AllocateSharedPtr(
        m_session, m_bndCondExpansions, m_bndConditions, *m_exp, m_graph,
        m_comm, true, "DefaultVar", ImpType);
 
    PeriodicMap periodicTraces = (m_expType == e1D)   ? m_periodicVerts
                                 : (m_expType == e2D) ? m_periodicEdges
                                                      : m_periodicFaces;
 
    m_traceMap = MemoryManager<AssemblyMapDG>::AllocateSharedPtr(
        m_session, m_graph, m_trace, *this, m_bndCondExpansions,
        m_bndConditions, periodicTraces, variable);
 
    if (m_session->DefinesCmdLineArgument("verbose"))
    {
        m_traceMap->PrintStats(std::cout, variable);
    }
 
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    // Scatter trace to elements. For each element, we find
    // the trace  associated to each elemental trace. The element
    // then retains a pointer to the elemental trace, to
    // ensure uniqueness of normals when retrieving from two
    // adjoining elements which do not lie in a plane.
    for (int i = 0; i < m_exp->size(); ++i)
    {
        for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
        {
            LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];
            ;
 
            exp->SetAdjacentElementExp(j, (*m_exp)[i]);
            (*m_exp)[i]->SetTraceExp(j, exp);
        }
    }
 
    // Set up physical normals
    SetUpPhysNormals();
 
    // Create interface exchange object
    m_interfaceMap =
        MemoryManager<InterfaceMapDG>::AllocateSharedPtr(m_graph, m_trace);
 
    int cnt, n;
 
    // Identify boundary trace
    for (cnt = 0, n = 0; n < m_bndCondExpansions.size(); ++n)
    {
        if (m_bndConditions[n]->GetBoundaryConditionType() !=
            SpatialDomains::ePeriodic)
        {
            for (int v = 0; v < m_bndCondExpansions[n]->GetExpSize(); ++v)
            {
                m_boundaryTraces.insert(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
            }
            cnt += m_bndCondExpansions[n]->GetExpSize();
        }
    }
 
    // Set up information for periodic boundary conditions.
    std::unordered_map<int, pair<int, int>> perTraceToExpMap;
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            auto it = periodicTraces.find((*m_exp)[n]->GetGeom()->GetTid(v));
 
            if (it != periodicTraces.end())
            {
                perTraceToExpMap[it->first] = make_pair(n, v);
            }
        }
    }
 
    // Set up left-adjacent tracelist.
    m_leftAdjacentTraces.resize(cnt);
 
    // count size of trace
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            m_leftAdjacentTraces[cnt] = IsLeftAdjacentTrace(n, v);
        }
    }
 
    // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
 
            // Check to see if this trace is periodic.
            auto it = periodicTraces.find(GeomId);
 
            if (it != periodicTraces.end())
            {
                const PeriodicEntity &ent = it->second[0];
                auto it2                  = perTraceToExpMap.find(ent.id);
 
                if (it2 == perTraceToExpMap.end())
                {
                    if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
                        !ent.isLocal)
                    {
                        continue;
                    }
                    else
                    {
                        ASSERTL1(false, "Periodic trace not found!");
                    }
                }
 
                ASSERTL1(m_leftAdjacentTraces[cnt],
                         "Periodic trace in non-forward space?");
 
                int offset = m_trace->GetPhys_Offset(
                    (m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
 
                int offset2 = m_trace->GetPhys_Offset(
                    (m_traceMap->GetElmtToTrace())[it2->second.first]
                                                  [it2->second.second]
                                                      ->GetElmtId());
 
                switch (m_expType)
                {
                    case e1D:
                    {
                        m_periodicFwdCopy.push_back(offset);
                        m_periodicBwdCopy.push_back(offset2);
                    }
                    break;
                    case e2D:
                    {
                        // Calculate relative orientations between trace to
                        // calculate copying map.
                        int nquad = elmtToTrace[n][v]->GetNumPoints(0);
 
                        vector<int> tmpBwd(nquad);
                        vector<int> tmpFwd(nquad);
 
                        if (ent.orient == StdRegions::eForwards)
                        {
                            for (int i = 0; i < nquad; ++i)
                            {
                                tmpBwd[i] = offset2 + i;
                                tmpFwd[i] = offset + i;
                            }
                        }
                        else
                        {
                            for (int i = 0; i < nquad; ++i)
                            {
                                tmpBwd[i] = offset2 + i;
                                tmpFwd[i] = offset + nquad - i - 1;
                            }
                        }
 
                        for (int i = 0; i < nquad; ++i)
                        {
                            m_periodicFwdCopy.push_back(tmpFwd[i]);
                            m_periodicBwdCopy.push_back(tmpBwd[i]);
                        }
                    }
                    break;
                    case e3D:
                    {
                        // Calculate relative orientations between faces to
                        // calculate copying map.
                        int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
                        int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
 
                        vector<int> tmpBwd(nquad1 * nquad2);
                        vector<int> tmpFwd(nquad1 * nquad2);
 
                        if (ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                            ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1; ++j)
                                {
                                    tmpBwd[i * nquad1 + j] =
                                        offset2 + i * nquad1 + j;
                                    tmpFwd[i * nquad1 + j] =
                                        offset + j * nquad2 + i;
                                }
                            }
                        }
                        else
                        {
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1; ++j)
                                {
                                    tmpBwd[i * nquad1 + j] =
                                        offset2 + i * nquad1 + j;
                                    tmpFwd[i * nquad1 + j] =
                                        offset + i * nquad1 + j;
                                }
                            }
                        }
 
                        if (ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            // Reverse x direction
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1 / 2; ++j)
                                {
                                    swap(tmpFwd[i * nquad1 + j],
                                         tmpFwd[i * nquad1 + nquad1 - j - 1]);
                                }
                            }
                        }
 
                        if (ent.orient ==
                                StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            // Reverse y direction
                            for (int j = 0; j < nquad1; ++j)
                            {
                                for (int i = 0; i < nquad2 / 2; ++i)
                                {
                                    swap(tmpFwd[i * nquad1 + j],
                                         tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
                                }
                            }
                        }
 
                        for (int i = 0; i < nquad1 * nquad2; ++i)
                        {
                            m_periodicFwdCopy.push_back(tmpFwd[i]);
                            m_periodicBwdCopy.push_back(tmpBwd[i]);
                        }
                    }
                    break;
                    default:
                        ASSERTL1(false, "not set up");
                }
            }
        }
    }
 
    m_locTraceToTraceMap = MemoryManager<LocTraceToTraceMap>::AllocateSharedPtr(
        *this, m_trace, elmtToTrace, m_leftAdjacentTraces);
}

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(), DisContField(), 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

Fwd	The trace quantities associated with left (fwd) adjancent elmt.
Bwd	The trace quantities associated with right (bwd) adjacent elet.
outarray	Resulting coefficient space.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3722 of file DisContField.cpp.

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

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

Fn	The trace quantities.
outarray	Resulting 3D coefficient space.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3551 of file DisContField.cpp.

{
    if (m_expType == e1D)
    {
        int n, offset, t_offset;
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace               = m_traceMap->GetElmtToTrace();
        vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
        for (n = 0; n < GetExpSize(); ++n)
        {
            // Number of coefficients on each element
            int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
            offset        = GetCoeff_Offset(n);
            // Implementation for every points except Gauss points
            if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
                LibUtilities::eGauss_Lagrange)
            {
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
                if (negatedFluxNormal[2 * n])
                {
                    outarray[offset] -= Fn[t_offset];
                }
                else
                {
                    outarray[offset] += Fn[t_offset];
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
                if (negatedFluxNormal[2 * n + 1])
                {
                    outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
                        Fn[t_offset];
                }
                else
                {
                    outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
                        Fn[t_offset];
                }
            }
            else
            {
                int j;
                static DNekMatSharedPtr m_Ixm, m_Ixp;
                static int sav_ncoeffs = 0;
                if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
                {
                    LibUtilities::BasisSharedPtr BASE;
                    const LibUtilities::PointsKey BS_p(
                        e_ncoeffs, LibUtilities::eGaussGaussLegendre);
                    const LibUtilities::BasisKey BS_k(
                        LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
                    BASE = LibUtilities::BasisManager()[BS_k];
                    Array<OneD, NekDouble> coords(1, 0.0);
                    coords[0]   = -1.0;
                    m_Ixm       = BASE->GetI(coords);
                    coords[0]   = 1.0;
                    m_Ixp       = BASE->GetI(coords);
                    sav_ncoeffs = e_ncoeffs;
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
                if (negatedFluxNormal[2 * n])
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] -=
                            (m_Ixm->GetPtr())[j] * Fn[t_offset];
                    }
                }
                else
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] +=
                            (m_Ixm->GetPtr())[j] * Fn[t_offset];
                    }
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
                if (negatedFluxNormal[2 * n + 1])
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] -=
                            (m_Ixp->GetPtr())[j] * Fn[t_offset];
                    }
                }
                else
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] +=
                            (m_Ixp->GetPtr())[j] * Fn[t_offset];
                    }
                }
            }
        }
    }
    else // other dimensions
    {
        // Basis definition on each element
        LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
        if ((m_expType != e1D) &&
            (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
        {
            Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
            m_trace->IProductWRTBase(Fn, Fcoeffs);
            m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
                                                              outarray);
        }
        else
        {
            int e, n, offset, t_offset;
            Array<OneD, NekDouble> e_outarray;
            Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
                &elmtToTrace = m_traceMap->GetElmtToTrace();
            if (m_expType == e2D)
            {
                for (n = 0; n < GetExpSize(); ++n)
                {
                    offset = GetCoeff_Offset(n);
                    for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
                    {
                        t_offset = GetTrace()->GetPhys_Offset(
                            elmtToTrace[n][e]->GetElmtId());
                        (*m_exp)[n]->AddEdgeNormBoundaryInt(
                            e, elmtToTrace[n][e], Fn + t_offset,
                            e_outarray = outarray + offset);
                    }
                }
            }
            else if (m_expType == e3D)
            {
                for (n = 0; n < GetExpSize(); ++n)
                {
                    offset = GetCoeff_Offset(n);
                    for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
                    {
                        t_offset = GetTrace()->GetPhys_Offset(
                            elmtToTrace[n][e]->GetElmtId());
                        (*m_exp)[n]->AddFaceNormBoundaryInt(
                            e, elmtToTrace[n][e], Fn + t_offset,
                            e_outarray = outarray + offset);
                    }
                }
            }
        }
    }
}

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

{
    Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
 
    m_trace->IProductWRTBase(FwdFlux, FCoeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
    m_trace->IProductWRTBase(BwdFlux, FCoeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
}

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

{
    // Basis definition on each element
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
    {
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
 
        Array<OneD, NekDouble> invals =
            tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
 
        m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
 
        m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
 
        m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
    }
    else
    {
        NEKERROR(ErrorUtil::efatal,
                 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
    }
}

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

v_EvaluateBoundaryConditions()

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

overrideprotectedvirtual

Evaluate all boundary conditions at a given time..

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3903 of file DisContField.cpp.

{
    int i;
    int npoints;
 
    MultiRegions::ExpListSharedPtr locExpList;
 
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        if (m_bndCondExpansions[i] &&
            (time == 0.0 || m_bndConditions[i]->IsTimeDependent()))
        {
            m_bndCondBndWeight[i] = 1.0;
            locExpList            = m_bndCondExpansions[i];
 
            npoints = locExpList->GetNpoints();
            Array<OneD, NekDouble> x0(npoints, 0.0);
            Array<OneD, NekDouble> x1(npoints, 0.0);
            Array<OneD, NekDouble> x2(npoints, 0.0);
 
            locExpList->GetCoords(x0, x1, x2);
 
            if (x2_in != NekConstants::kNekUnsetDouble &&
                x3_in != NekConstants::kNekUnsetDouble)
            {
                Vmath::Fill(npoints, x2_in, x1, 1);
                Vmath::Fill(npoints, x3_in, x2, 1);
            }
            else if (x2_in != NekConstants::kNekUnsetDouble)
            {
                Vmath::Fill(npoints, x2_in, x2, 1);
            }
 
            // treat 1D expansions separately since we only
            // require an evaluation at a point rather than
            // any projections or inner products that are not
            // available in a PointExp
            if (m_expType == e1D)
            {
                if (m_bndConditions[i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::DirichletBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_dirichletCondition)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                        m_bndCondExpansions[i]->SetPhys(
                            n, m_bndCondExpansions[i]->GetCoeff(n));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::NeumannBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_neumannCondition)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::RobinBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_robinFunction)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::ePeriodic)
                {
                    continue;
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNotDefined)
                {
                }
                else
                {
                    NEKERROR(ErrorUtil::efatal,
                             "This type of BC not implemented yet");
                }
            }
            else // 2D and 3D versions
            {
                if (m_bndConditions[i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    SpatialDomains::DirichletBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::DirichletBoundaryCondition>(
                            m_bndConditions[i]);
 
                    Array<OneD, NekDouble> valuesFile(npoints, 1.0),
                        valuesExp(npoints, 1.0);
 
                    string bcfilename = bcPtr->m_filename;
                    string exprbcs    = bcPtr->m_expr;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                        valuesFile = locExpList->GetPhys();
                    }
 
                    if (exprbcs != "")
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::DirichletBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_dirichletCondition;
 
                        condition->Evaluate(x0, x1, x2, time, valuesExp);
                    }
 
                    Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
                                locExpList->UpdatePhys(), 1);
                    locExpList->FwdTransBndConstrained(
                        locExpList->GetPhys(), locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann)
                {
                    SpatialDomains::NeumannBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::NeumannBoundaryCondition>(
                            m_bndConditions[i]);
                    string bcfilename = bcPtr->m_filename;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                    }
                    else
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::NeumannBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_neumannCondition;
                        condition->Evaluate(x0, x1, x2, time,
                                            locExpList->UpdatePhys());
                    }
 
                    locExpList->IProductWRTBase(locExpList->GetPhys(),
                                                locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
                {
                    SpatialDomains::RobinBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::RobinBoundaryCondition>(
                            m_bndConditions[i]);
 
                    string bcfilename = bcPtr->m_filename;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                    }
                    else
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::RobinBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_robinFunction;
                        condition->Evaluate(x0, x1, x2, time,
                                            locExpList->UpdatePhys());
                    }
 
                    locExpList->IProductWRTBase(locExpList->GetPhys(),
                                                locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::ePeriodic)
                {
                    continue;
                }
                else
                {
                    NEKERROR(ErrorUtil::efatal,
                             "This type of BC not implemented yet");
                }
            }
        }
    }
}

References Nektar::MultiRegions::e1D, Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eNotDefined, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, 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 3470 of file DisContField.cpp.

{
    ASSERTL1(m_physState == true, "local physical space is not true ");
    v_ExtractTracePhys(m_phys, outarray);
}

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, bool  gridVelocity = false  )

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

inarray	An array containing the 1D data from which we wish to extract the edge data.
outarray	The resulting edge information.
gridVelocity	Avoid performing parallel exchanges of the grid velocity

This will not work for non-boundary expansions

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3491 of file DisContField.cpp.

{
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
    {
        Vmath::Zero(outarray.size(), outarray, 1);
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNFwdLocTracePts());
        m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
        m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
        if (!gridVelocity)
        {
            m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray,
                                                             outarray);
        }
    }
    else
    {
        // Loop over elemente and collect forward expansion
        int nexp = GetExpSize();
        int n, p, offset, phys_offset;
        Array<OneD, NekDouble> t_tmp;
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
        ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
                 "input array is of insufficient length");
        for (n = 0; n < nexp; ++n)
        {
            phys_offset = GetPhys_Offset(n);
            for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
            {
                offset =
                    m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
                (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
                                              inarray + phys_offset,
                                              t_tmp = outarray + offset);
            }
        }
    }
}

References ASSERTL1, Nektar::LibUtilities::eGauss_Lagrange, Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_locTraceToTraceMap, m_trace, m_traceMap, 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 3314 of file DisContField.cpp.

{
    FillBwdWithBoundCond(Fwd, Bwd, m_bndConditions, m_bndCondExpansions,
                         PutFwdInBwdOnBCs);
}

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

{
    int cnt;
    int npts;
    int e = 0;
 
    // Fill boundary conditions into missing elements
    int id2 = 0;
 
    for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
    {
 
        if (m_bndConditions[n]->GetBoundaryConditionType() ==
            SpatialDomains::eDirichlet)
        {
            for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
            {
                npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                id2  = m_trace->GetPhys_Offset(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
                Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
            }
 
            cnt += e;
        }
        else if (m_bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann ||
                 m_bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
        {
            for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
            {
                npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                id2  = m_trace->GetPhys_Offset(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
                Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
            }
 
            cnt += e;
        }
        else if (m_bndConditions[n]->GetBoundaryConditionType() !=
                 SpatialDomains::ePeriodic)
        {
            NEKERROR(ErrorUtil::efatal,
                     "Method not set up for this boundary condition.");
        }
    }
}

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

{
    return m_bndCondBndWeight;
}

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

{
    return m_bndCondExpansions;
}

References m_bndCondExpansions.

v_GetBndConditions()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3025 of file DisContField.cpp.

{
    return m_bndConditions;
}

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

{
    int n, cnt;
    std::vector<unsigned int> eIDs;
 
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    // Skip other boundary regions
    for (cnt = n = 0; n < i; ++n)
    {
        cnt += m_bndCondExpansions[n]->GetExpSize();
    }
 
    // Populate eIDs with information from BoundaryToElmtMap
    for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
    {
        eIDs.push_back(ElmtID[cnt + n]);
    }
 
    // Create expansion list
    result = MemoryManager<ExpList>::AllocateSharedPtr(
        *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
    // Copy phys and coeffs to new explist
    if (DeclareCoeffPhysArrays)
    {
        int nq;
        int offsetOld, offsetNew;
        Array<OneD, NekDouble> tmp1, tmp2;
        for (n = 0; n < result->GetExpSize(); ++n)
        {
            nq        = GetExp(ElmtID[cnt + n])->GetTotPoints();
            offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
            offsetNew = result->GetPhys_Offset(n);
            Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
                         tmp2 = result->UpdatePhys() + offsetNew, 1);
            nq        = GetExp(ElmtID[cnt + n])->GetNcoeffs();
            offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
            offsetNew = result->GetCoeff_Offset(n);
            Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
                         tmp2 = result->UpdateCoeffs() + offsetNew, 1);
        }
    }
}

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

{
 
    if (m_BCtoElmMap.size() == 0)
    {
        switch (m_expType)
        {
            case e1D:
            {
                map<int, int> VertGID;
                int i, n, id;
                int bid, cnt, Vid;
                int nbcs = 0;
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    nbcs += m_bndCondExpansions[i]->GetExpSize();
                }
 
                // make sure arrays are of sufficient length
                m_BCtoElmMap   = Array<OneD, int>(nbcs, -1);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                // setup map of all global ids along boundary
                cnt = 0;
                for (n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i)
                    {
                        Vid = m_bndCondExpansions[n]
                                  ->GetExp(i)
                                  ->GetGeom()
                                  ->GetVid(0);
                        VertGID[Vid] = cnt++;
                    }
                }
 
                // Loop over elements and find verts that match;
                LocalRegions::ExpansionSharedPtr exp;
                for (cnt = n = 0; n < GetExpSize(); ++n)
                {
                    exp = (*m_exp)[n];
                    for (i = 0; i < exp->GetNverts(); ++i)
                    {
                        id = exp->GetGeom()->GetVid(i);
 
                        if (VertGID.count(id) > 0)
                        {
                            bid = VertGID.find(id)->second;
                            ASSERTL1(m_BCtoElmMap[bid] == -1,
                                     "Edge already set");
                            m_BCtoElmMap[bid]   = n;
                            m_BCtoTraceMap[bid] = i;
                            cnt++;
                        }
                    }
                }
                ASSERTL1(cnt == nbcs,
                         "Failed to visit all boundary condtiions");
            }
            break;
            case e2D:
            {
                map<int, int> globalIdMap;
                int i, n;
                int cnt;
                int nbcs = 0;
 
                // Populate global ID map (takes global geometry
                // ID to local expansion list ID).
                for (i = 0; i < GetExpSize(); ++i)
                {
                    globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                }
 
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    nbcs += m_bndCondExpansions[i]->GetExpSize();
                }
 
                // Initialize arrays
                m_BCtoElmMap   = Array<OneD, int>(nbcs);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                LocalRegions::Expansion1DSharedPtr exp1d;
                cnt = 0;
                for (n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
                         ++i, ++cnt)
                    {
                        exp1d = m_bndCondExpansions[n]
                                    ->GetExp(i)
                                    ->as<LocalRegions::Expansion1D>();
 
                        // Use edge to element map from MeshGraph.
                        SpatialDomains::GeometryLinkSharedPtr tmp =
                            m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
 
                        m_BCtoElmMap[cnt] =
                            globalIdMap[(*tmp)[0].first->GetGlobalID()];
                        m_BCtoTraceMap[cnt] = (*tmp)[0].second;
                    }
                }
            }
            break;
            case e3D:
            {
                map<int, int> globalIdMap;
                int i, n;
                int cnt;
                int nbcs = 0;
 
                // Populate global ID map (takes global geometry ID to local
                // expansion list ID).
                LocalRegions::Expansion3DSharedPtr exp3d;
                for (i = 0; i < GetExpSize(); ++i)
                {
                    globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                }
 
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    if (m_bndCondExpansions[i])
                    {
                        nbcs += m_bndCondExpansions[i]->GetExpSize();
                    }
                }
 
                // Initialize arrays
                m_BCtoElmMap   = Array<OneD, int>(nbcs);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                LocalRegions::Expansion2DSharedPtr exp2d;
                for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
                         ++i, ++cnt)
                    {
                        exp2d = m_bndCondExpansions[n]
                                    ->GetExp(i)
                                    ->as<LocalRegions::Expansion2D>();
 
                        SpatialDomains::GeometryLinkSharedPtr tmp =
                            m_graph->GetElementsFromFace(exp2d->GetGeom2D());
                        m_BCtoElmMap[cnt] =
                            globalIdMap[tmp->at(0).first->GetGlobalID()];
                        m_BCtoTraceMap[cnt] = tmp->at(0).second;
                    }
                }
            }
            break;
            default:
                ASSERTL1(false, "Not setup for this expansion");
                break;
        }
    }
 
    ElmtID  = m_BCtoElmMap;
    TraceID = m_BCtoTraceMap;
}

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

{
    v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
}

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

{
    // Is this zeroing necessary?
    // Zero forward/backward vectors.
    Vmath::Zero(Fwd.size(), Fwd, 1);
    Vmath::Zero(Bwd.size(), Bwd, 1);
 
    // Basis definition on each element
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
    {
        // blocked routine
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNLocTracePts());
 
        m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
        m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
 
        Array<OneD, NekDouble> invals =
            tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
        m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
    }
    else
    {
        // Loop over elements and collect forward expansion
        auto nexp = GetExpSize();
        Array<OneD, NekDouble> e_tmp;
        LocalRegions::ExpansionSharedPtr exp;
 
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
 
        for (int n = 0, cnt = 0; n < nexp; ++n)
        {
            exp              = (*m_exp)[n];
            auto phys_offset = GetPhys_Offset(n);
 
            for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
            {
                auto offset =
                    m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
 
                e_tmp =
                    (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
 
                exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
                                      e_tmp);
            }
        }
    }
 
    DisContField::v_PeriodicBwdCopy(Fwd, Bwd);
 
    if (FillBnd)
    {
        v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
    }
 
    if (DoExchange)
    {
        // Do parallel exchange for forwards/backwards spaces.
        m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
 
        // Do exchange of interface traces (local and parallel)
        // We may have to split this out into separate local and
        // parallel for IP method???
        m_interfaceMap->ExchangeTrace(Fwd, Bwd);
    }
}

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

v_GetInterfaceMap()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3000 of file DisContField.cpp.

{
    return m_interfaceMap;
}

References m_interfaceMap.

v_GetLeftAdjacentTraces()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3014 of file DisContField.cpp.

{
    return m_leftAdjacentTraces;
}

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

{
    if (locTraceFwd != NullNekDouble1DArray)
    {
        m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
                                                              locTraceFwd);
    }
 
    if (locTraceBwd != NullNekDouble1DArray)
    {
        m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
                                                              locTraceBwd);
    }
}

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

{
    return m_locTraceToTraceMap;
}

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

{
    periodicVerts = m_periodicVerts;
    periodicEdges = m_periodicEdges;
    periodicFaces = m_periodicFaces;
}

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

{
    int i, cnt;
    map<int, RobinBCInfoSharedPtr> returnval;
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        MultiRegions::ExpListSharedPtr locExpList;
 
        if (m_bndConditions[i]->GetBoundaryConditionType() ==
            SpatialDomains::eRobin)
        {
            int e, elmtid;
            Array<OneD, NekDouble> Array_tmp;
 
            locExpList = m_bndCondExpansions[i];
 
            int npoints = locExpList->GetNpoints();
            Array<OneD, NekDouble> x0(npoints, 0.0);
            Array<OneD, NekDouble> x1(npoints, 0.0);
            Array<OneD, NekDouble> x2(npoints, 0.0);
            Array<OneD, NekDouble> coeffphys(npoints);
 
            locExpList->GetCoords(x0, x1, x2);
 
            LibUtilities::EquationSharedPtr coeffeqn =
                std::static_pointer_cast<
                    SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i])
                    ->m_robinPrimitiveCoeff;
 
            // evalaute coefficient
            coeffeqn->Evaluate(x0, x1, x2, 0.0, coeffphys);
 
            for (e = 0; e < locExpList->GetExpSize(); ++e)
            {
                RobinBCInfoSharedPtr rInfo =
                    MemoryManager<RobinBCInfo>::AllocateSharedPtr(
                        TraceID[cnt + e],
                        Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
 
                elmtid = ElmtID[cnt + e];
                // make link list if necessary
                if (returnval.count(elmtid) != 0)
                {
                    rInfo->next = returnval.find(elmtid)->second;
                }
                returnval[elmtid] = rInfo;
            }
        }
        cnt += m_bndCondExpansions[i]->GetExpSize();
    }
 
    return returnval;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SpatialDomains::eRobin, 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 2985 of file DisContField.cpp.

{
    if (m_trace == NullExpListSharedPtr)
    {
        SetUpDG();
    }
 
    return m_trace;
}

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

{
    return m_traceMap;
}

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 StdRegions::VarFactorsMap &  varfactors, const Array< OneD, const NekDouble > &  dirForcing, const bool  PhysSpaceForcing  )

overrideprotectedvirtual

Solve the Helmholtz equation.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3735 of file DisContField.cpp.

{
    int i, n, cnt, nbndry;
    int nexp = GetExpSize();
    Array<OneD, NekDouble> f(m_ncoeffs);
    DNekVec F(m_ncoeffs, f, eWrapper);
    Array<OneD, NekDouble> e_f, e_l;
    //----------------------------------
    // Setup RHS Inner product if required
    //----------------------------------
    if (PhysSpaceForcing)
    {
        IProductWRTBase(inarray, f);
        Vmath::Neg(m_ncoeffs, f, 1);
    }
    else
    {
        Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
    }
    //----------------------------------
    // Solve continuous Boundary System
    //----------------------------------
    int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
    int NumDirBCs  = m_traceMap->GetNumLocalDirBndCoeffs();
    int e_ncoeffs;
    GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
                                 NullAssemblyMapSharedPtr, factors, varcoeff);
    const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
    // Retrieve number of local trace space coefficients N_{\lambda},
    // and set up local elemental trace solution \lambda^e.
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
    Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
    DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
    //----------------------------------
    // Evaluate Trace Forcing
    // Kirby et al, 2010, P23, Step 5.
    //----------------------------------
    // Determing <u_lam,f> terms using HDGLamToU matrix
    for (cnt = n = 0; n < nexp; ++n)
    {
        nbndry    = (*m_exp)[n]->NumDGBndryCoeffs();
        e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
        e_f       = f + m_coeff_offset[n];
        e_l       = bndrhs + cnt;
        // use outarray as tmp space
        DNekVec Floc(nbndry, e_l, eWrapper);
        DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
        Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
        cnt += nbndry;
    }
    Array<OneD, const int> bndCondMap =
        m_traceMap->GetBndCondCoeffsToLocalTraceMap();
    Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
    // Copy Dirichlet boundary conditions and weak forcing
    // into trace space
    int locid;
    cnt = 0;
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        const Array<OneD, const NekDouble> bndcoeffs =
            m_bndCondExpansions[i]->GetCoeffs();
 
        if (m_bndConditions[i]->GetBoundaryConditionType() ==
            SpatialDomains::eDirichlet)
        {
            for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
            {
                locid            = bndCondMap[cnt + j];
                loclambda[locid] = Sign[locid] * bndcoeffs[j];
            }
        }
        else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann ||
                 m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
        {
            // Add weak boundary condition to trace forcing
            for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
            {
                locid = bndCondMap[cnt + j];
                bndrhs[locid] += Sign[locid] * bndcoeffs[j];
            }
        }
        else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                 SpatialDomains::ePeriodic)
        {
            // skip increment of cnt if ePeriodic
            // because bndCondMap does not include ePeriodic
            continue;
        }
        cnt += (m_bndCondExpansions[i])->GetNcoeffs();
    }
 
    //----------------------------------
    // Solve trace problem: \Lambda = K^{-1} F
    // K is the HybridDGHelmBndLam matrix.
    //----------------------------------
    if (GloBndDofs - NumDirBCs > 0)
    {
        GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam, m_traceMap,
                            factors, varcoeff);
        GlobalLinSysSharedPtr LinSys = GetGlobalBndLinSys(key);
        LinSys->Solve(bndrhs, loclambda, m_traceMap);
 
        // For consistency with previous version put global
        // solution into m_trace->m_coeffs
        m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
    }
 
    //----------------------------------
    // Internal element solves
    //----------------------------------
    GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
                                  NullAssemblyMapSharedPtr, factors, varcoeff);
 
    const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
    DNekVec out(m_ncoeffs, outarray, eWrapper);
    Vmath::Zero(m_ncoeffs, outarray, 1);
 
    //  out =  u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
    out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
 
    // Return empty GlobalLinSysKey
    return NullGlobalLinSysKey;
}

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

v_PeriodicBwdCopy()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3042 of file DisContField.cpp.

{
    for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
    {
        Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
    }
}

References m_periodicBwdCopy, and m_periodicFwdCopy.

Referenced by GetFwdBwdTracePhys(), and v_GetFwdBwdTracePhys().

v_PeriodicBwdRot()

void Nektar::MultiRegions::DisContField::v_PeriodicBwdRot ( Array< OneD, Array< OneD, NekDouble > > &  Bwd )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3051 of file DisContField.cpp.

{
    int cnt = 0;
    for (int n = 0; n < m_bndConditions.size(); ++n)
    {
        // check to see if boundary is rotationally aligned
        if (boost::icontains(m_bndConditions[n]->GetUserDefined(), "Rotated"))
        {
            vector<string> tmpstr;
 
            boost::split(tmpstr, m_bndConditions[n]->GetUserDefined(),
                         boost::is_any_of(":"));
 
            ASSERTL1(tmpstr.size() > 2,
                     "Expected Rotated user defined string to "
                     "contain direction and rotation angle "
                     "and optionally a tolerance, "
                     "i.e. Rotated:dir:PI/2:1e-6");
 
            ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                         (tmpstr[1] == "z"),
                     "Rotated Dir is "
                     "not specified as x,y or z");
 
            RotPeriodicInfo RotInfo;
            RotInfo.m_dir = (tmpstr[1] == "x") ? 0 : (tmpstr[1] == "y") ? 1 : 2;
 
            LibUtilities::Interpreter strEval;
            int ExprId      = strEval.DefineFunction("", tmpstr[2]);
            RotInfo.m_angle = strEval.Evaluate(ExprId);
 
            auto ne = m_bndCondExpansions[n]->GetExpSize();
 
            // Loop over each element in the boundary condition and rotate
            // velocity
            for (int e = 0; e < ne; ++e)
            {
                auto id2 = m_trace->GetPhys_Offset(
                    m_traceMap->GetPerBndCondIDToGlobalTraceID(e + cnt));
                int npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                Rotate(Bwd, RotInfo.m_dir, -RotInfo.m_angle, id2, npts);
            }
            cnt += ne;
        }
    }
}

References ASSERTL1, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::LibUtilities::Interpreter::Evaluate(), Nektar::MultiRegions::RotPeriodicInfo::m_angle, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::RotPeriodicInfo::m_dir, m_trace, m_traceMap, and Rotate().

v_PeriodicDeriveBwdRot()

void Nektar::MultiRegions::DisContField::v_PeriodicDeriveBwdRot ( TensorOfArray3D< NekDouble > &  Bwd )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3098 of file DisContField.cpp.

{
    int cnt = 0;
    for (int n = 0; n < m_bndConditions.size(); ++n)
    {
        // check to see if boundary is rotationally aligned
        if (boost::icontains(m_bndConditions[n]->GetUserDefined(), "Rotated"))
        {
            vector<string> tmpstr;
 
            boost::split(tmpstr, m_bndConditions[n]->GetUserDefined(),
                         boost::is_any_of(":"));
 
            ASSERTL1(tmpstr.size() > 2,
                     "Expected Rotated user defined string to "
                     "contain direction and rotation angle "
                     "and optionally a tolerance, "
                     "i.e. Rotated:dir:PI/2:1e-6");
 
            ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                         (tmpstr[1] == "z"),
                     "Rotated Dir is "
                     "not specified as x,y or z");
 
            RotPeriodicInfo RotInfo;
            RotInfo.m_dir = (tmpstr[1] == "x") ? 0 : (tmpstr[1] == "y") ? 1 : 2;
 
            LibUtilities::Interpreter strEval;
            int ExprId      = strEval.DefineFunction("", tmpstr[2]);
            RotInfo.m_angle = strEval.Evaluate(ExprId);
 
            auto ne = m_bndCondExpansions[n]->GetExpSize();
            // Loop over each element in the boundary condition and rotate
            // velocity
            for (int e = 0; e < ne; ++e)
            {
                auto id2 = m_trace->GetPhys_Offset(
                    m_traceMap->GetPerBndCondIDToGlobalTraceID(e + cnt));
                int npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                DeriveRotate(Bwd, RotInfo.m_dir, -RotInfo.m_angle, id2, npts);
            }
            cnt += ne;
        }
    }
}

References ASSERTL1, Nektar::LibUtilities::Interpreter::DefineFunction(), DeriveRotate(), Nektar::LibUtilities::Interpreter::Evaluate(), Nektar::MultiRegions::RotPeriodicInfo::m_angle, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::RotPeriodicInfo::m_dir, m_trace, and m_traceMap.

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

{
    ExpList::v_Reset();
 
    // Reset boundary condition expansions.
    for (int n = 0; n < m_bndCondExpansions.size(); ++n)
    {
        m_bndCondExpansions[n]->Reset();
    }
}

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

v_RotLocalBwdDeriveTrace()

void Nektar::MultiRegions::DisContField::v_RotLocalBwdDeriveTrace ( TensorOfArray3D< NekDouble > &  Bwd )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3154 of file DisContField.cpp.

{
    int nDim = GetCoordim(0);
    for (auto &interfaceTrace : m_interfaceMap->GetLocalInterface())
    {
 
        interfaceTrace->RotLocalBwdDeriveTrace(Bwd, nDim);
    }
}

References Nektar::MultiRegions::ExpList::GetCoordim(), and m_interfaceMap.

v_RotLocalBwdTrace()

void Nektar::MultiRegions::DisContField::v_RotLocalBwdTrace ( Array< OneD, Array< OneD, NekDouble > > &  Bwd )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3144 of file DisContField.cpp.

{
    int nDim = GetCoordim(0);
    for (auto &interfaceTrace : m_interfaceMap->GetLocalInterface())
    {
 
        interfaceTrace->RotLocalBwdTrace(Bwd, nDim);
    }
}

References Nektar::MultiRegions::ExpList::GetCoordim(), and m_interfaceMap.

v_SetBCsToHomogeneous()

void Nektar::MultiRegions::DisContField::v_SetBCsToHomogeneous ( void  )

overrideprotectedvirtual

Set boundary conditions to be homogeneous.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4121 of file DisContField.cpp.

{
    int i;
    int npoints;
 
    MultiRegions::ExpListSharedPtr locExpList;
 
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        locExpList = m_bndCondExpansions[i];
 
        npoints = locExpList->GetNpoints();
 
        // treat 1D expansions separately since we only
        // require an evaluation at a point rather than
        // any projections or inner products that are not
        // available in a PointExp
        if (m_expType == e1D)
        {
            NEKERROR(ErrorUtil::efatal,
                     "Function needs setting up for 1D boundary expansion");
        }
        else // 2D and 3D versions
        {
            if (m_bndConditions[i]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                SpatialDomains::DirichletBCShPtr bcPtr =
                    std::static_pointer_cast<
                        SpatialDomains::DirichletBoundaryCondition>(
                        m_bndConditions[i]);
 
                Array<OneD, NekDouble> valuesExp(npoints, 0.0);
                bcPtr->m_filename = "";
                bcPtr->m_expr     = "0.0";
 
                bcPtr->m_dirichletCondition =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), bcPtr->m_expr);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann)
            {
                SpatialDomains::NeumannBCShPtr bcPtr = std::static_pointer_cast<
                    SpatialDomains::NeumannBoundaryCondition>(
                    m_bndConditions[i]);
 
                bcPtr->m_filename = "";
                std::string eqn   = "0.0";
 
                bcPtr->m_neumannCondition =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), eqn);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
            {
                SpatialDomains::RobinBCShPtr bcPtr = std::static_pointer_cast<
                    SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i]);
 
                bcPtr->m_filename = "";
                std::string eqn   = "0.0";
 
                bcPtr->m_robinFunction =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), eqn);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::ePeriodic)
            {
                continue;
            }
            else
            {
                NEKERROR(ErrorUtil::efatal,
                         "This type of BC not implemented yet");
            }
        }
    }
}

References Nektar::MultiRegions::e1D, Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_expType, Nektar::SpatialDomains::NeumannBoundaryCondition::m_filename, Nektar::SpatialDomains::RobinBoundaryCondition::m_filename, Nektar::MultiRegions::ExpList::m_session, NEKERROR, and Vmath::Zero().

v_SetBndCondBwdWeight()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3438 of file DisContField.cpp.

{
    m_bndCondBndWeight[index] = value;
}

References m_bndCondBndWeight.

v_UpdateBndCondExpansion()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3031 of file DisContField.cpp.

{
    return m_bndCondExpansions[i];
}

References m_bndCondExpansions.

v_UpdateBndConditions()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3036 of file DisContField.cpp.

{
    return m_bndConditions;
}

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

Definition at line 174 of file DisContField.h.

Referenced by GenerateBoundaryConditionExpansion(), GenerateBoundaryConditionExpansion(), v_EvaluateBoundaryConditions(), v_FillBwdWithBwdWeight(), v_GetBndCondBwdWeight(), and v_SetBndCondBwdWeight().

m_bndCondExpansions

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

protected

An object which contains the discretised boundary conditions.

It is an array of size equal to the number of boundary regions and consists of entries of the type MultiRegions::ExpList. Every entry corresponds to the spectral/hp expansion on a single boundary region. The values of the boundary conditions are stored as the coefficients of the one-dimensional expansion.

Definition at line 172 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), Nektar::MultiRegions::ContField::ContField(), DisContField(), DisContField(), FillBwdWithBoundCond(), GenerateBoundaryConditionExpansion(), GenerateBoundaryConditionExpansion(), Nektar::MultiRegions::ContField::LaplaceSolve(), SameTypeOfBoundaryConditions(), SetUpDG(), v_EvaluateBoundaryConditions(), Nektar::MultiRegions::ContField::v_FillBndCondFromField(), Nektar::MultiRegions::ContField::v_FillBndCondFromField(), v_FillBwdWithBoundCond(), v_FillBwdWithBwdWeight(), Nektar::MultiRegions::ContField::v_GetBndCondExpansions(), v_GetBndCondExpansions(), v_GetBndElmtExpansion(), v_GetBoundaryToElmtMap(), v_GetRobinBCInfo(), Nektar::MultiRegions::ContField::v_HelmSolve(), v_HelmSolve(), Nektar::MultiRegions::ContField::v_ImposeDirichletConditions(), Nektar::MultiRegions::ContField::v_ImposeNeumannConditions(), Nektar::MultiRegions::ContField::v_ImposeRobinConditions(), Nektar::MultiRegions::ContField::v_LinearAdvectionDiffusionReactionSolve(), Nektar::MultiRegions::ContField::v_LinearAdvectionReactionSolve(), v_PeriodicBwdRot(), v_PeriodicDeriveBwdRot(), v_Reset(), v_SetBCsToHomogeneous(), and v_UpdateBndCondExpansion().

m_bndConditions

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

protected

An array which contains the information about the boundary condition structure definition on the different boundary regions.

Definition at line 159 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), Nektar::MultiRegions::ContField::ContField(), GenerateBoundaryConditionExpansion(), GenerateBoundaryConditionExpansion(), Nektar::MultiRegions::ContField::LaplaceSolve(), SameTypeOfBoundaryConditions(), SetUpDG(), v_EvaluateBoundaryConditions(), v_FillBwdWithBoundCond(), v_FillBwdWithBwdWeight(), Nektar::MultiRegions::ContField::v_GetBndConditions(), v_GetBndConditions(), v_GetBoundaryToElmtMap(), v_GetRobinBCInfo(), Nektar::MultiRegions::ContField::v_HelmSolve(), v_HelmSolve(), Nektar::MultiRegions::ContField::v_ImposeDirichletConditions(), Nektar::MultiRegions::ContField::v_ImposeNeumannConditions(), Nektar::MultiRegions::ContField::v_ImposeRobinConditions(), Nektar::MultiRegions::ContField::v_LinearAdvectionDiffusionReactionSolve(), Nektar::MultiRegions::ContField::v_LinearAdvectionReactionSolve(), v_PeriodicBwdRot(), v_PeriodicDeriveBwdRot(), v_SetBCsToHomogeneous(), and v_UpdateBndConditions().

m_boundaryTraces

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

protected

A set storing the global IDs of any boundary Verts.

Definition at line 194 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 180 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 177 of file DisContField.h.

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

m_leftAdjacentTraces

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

protected

Definition at line 223 of file DisContField.h.

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

m_locElmtTrace

MultiRegions::ExpListSharedPtr Nektar::MultiRegions::DisContField::m_locElmtTrace

protected

Local Elemental trace expansions.

Definition at line 186 of file DisContField.h.

Referenced by DisContField(), DisContField(), and GetLocElmtTrace().

m_locTraceToTraceMap

LocTraceToTraceMapSharedPtr Nektar::MultiRegions::DisContField::m_locTraceToTraceMap

protected

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

Definition at line 229 of file DisContField.h.

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

m_negatedFluxNormal

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

private

Definition at line 381 of file DisContField.h.

Referenced by GetNegatedFluxNormal().

m_periodicBwdCopy

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

protected

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

Referenced by Nektar::MultiRegions::ContField::ContField(), 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 209 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), 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 216 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 199 of file DisContField.h.

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

m_trace

ExpListSharedPtr Nektar::MultiRegions::DisContField::m_trace

protected

Trace space storage for points between elements.

Definition at line 183 of file DisContField.h.

Referenced by DisContField(), DisContField(), EvaluateHDGPostProcessing(), FillBwdWithBoundCond(), L2_DGDeriv(), SetUpDG(), v_AddFwdBwdTraceIntegral(), v_AddTraceIntegral(), v_AddTraceIntegralToOffDiag(), v_ExtractTracePhys(), v_FillBwdWithBwdWeight(), v_GetFwdBwdTracePhys(), v_GetTrace(), v_HelmSolve(), v_PeriodicBwdRot(), and v_PeriodicDeriveBwdRot().

m_traceMap

AssemblyMapDGSharedPtr Nektar::MultiRegions::DisContField::m_traceMap

protected

Local to global DG mapping for trace space.

Definition at line 189 of file DisContField.h.

Referenced by DisContField(), EvaluateHDGPostProcessing(), FillBwdWithBoundCond(), GetFwdBwdTracePhys(), GetGlobalBndLinSys(), GetNegatedFluxNormal(), GetTraceElmtId(), IsLeftAdjacentTrace(), L2_DGDeriv(), SetUpDG(), v_AddTraceIntegral(), v_ExtractTracePhys(), v_FillBwdWithBwdWeight(), v_GetFwdBwdTracePhys(), v_GetTraceMap(), v_HelmSolve(), v_PeriodicBwdRot(), and v_PeriodicDeriveBwdRot().