Nektar++
Public Member Functions | List of all members
Nektar::CoupledLocalToGlobalC0ContMap Class Reference

#include <CoupledLocalToGlobalC0ContMap.h>

Inheritance diagram for Nektar::CoupledLocalToGlobalC0ContMap:
[legend]

Public Member Functions

 CoupledLocalToGlobalC0ContMap (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const MultiRegions::ExpListSharedPtr &pressure, const size_t nz_loc, const bool CheeckForSingularSys=true)
 
void FindEdgeIdToAddMeanPressure (std::vector< std::map< int, int > > &ReorderedGraphVertId, size_t &nel, const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId, size_t &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
 
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &comm, const std::string variable="DefaultVar")
 Default constructor. More...
 
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const BndCondExp &bndCondExp=NullExpListSharedPtrArray, const BndCond &bndConditions=SpatialDomains::NullBoundaryConditionShPtrArray, const bool checkIfSingular=false, const std::string variable="defaultVar", const PeriodicMap &periodicVerts=NullPeriodicMap, const PeriodicMap &periodicEdges=NullPeriodicMap, const PeriodicMap &periodicFaces=NullPeriodicMap)
 General constructor for expansions of all dimensions without boundary conditions. More...
 
 ~AssemblyMapCG () override
 Destructor. More...
 
std::set< ExtraDirDof > & GetCopyLocalDirDofs ()
 
std::set< int > & GetParallelDirBndSign ()
 
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
 AssemblyMap ()
 Default constructor. More...
 
 AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &comm, const std::string variable="DefaultVar")
 Constructor with a communicator. More...
 
 AssemblyMap (AssemblyMap *oldLevelMap, const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
 Constructor for next level in multi-level static condensation. More...
 
virtual ~AssemblyMap ()
 Destructor. More...
 
LibUtilities::CommSharedPtr GetComm ()
 Retrieves the communicator. More...
 
std::string GetVariable ()
 Retrieves the variable string. More...
 
size_t GetHash () const
 Retrieves the hash of this map. More...
 
int GetLocalToGlobalMap (const int i) const
 
int GetGlobalToUniversalMap (const int i) const
 
int GetGlobalToUniversalMapUnique (const int i) const
 
const Array< OneD, const int > & GetLocalToGlobalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMapUnique ()
 
NekDouble GetLocalToGlobalSign (const int i) const
 
const Array< OneD, NekDouble > & GetLocalToGlobalSign () const
 
void LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
 
void LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm=true) const
 
void GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
void GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
 
void Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
 
void UniversalAssemble (NekVector< NekDouble > &pGlobal) const
 
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
 
void UniversalAbsMaxBnd (Array< OneD, NekDouble > &bndvals)
 
int GetLocalToGlobalBndMap (const int i) const
 Retrieve the global index of a given local boundary mode. More...
 
const Array< OneD, const int > & GetLocalToGlobalBndMap ()
 Retrieve the global indices of the local boundary modes. More...
 
const Array< OneD, const int > & GetGlobalToUniversalBndMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalBndMapUnique ()
 
bool GetSignChange ()
 Returns true if using a modal expansion requiring a change of sign of some modes. More...
 
NekDouble GetLocalToGlobalBndSign (const int i) const
 Retrieve the sign change of a given local boundary mode. More...
 
Array< OneD, const NekDouble > GetLocalToGlobalBndSign () const
 Retrieve the sign change for all local boundary modes. More...
 
const Array< OneD, const int > & GetBndCondCoeffsToLocalCoeffsMap ()
 Retrieves the local indices corresponding to the boundary expansion modes. More...
 
const Array< OneD, NekDouble > & GetBndCondCoeffsToLocalCoeffsSign ()
 Returns the modal sign associated with a given boundary expansion mode. More...
 
const Array< OneD, const int > & GetBndCondCoeffsToLocalTraceMap ()
 Retrieves the local indices corresponding to the boundary expansion modes to global trace. More...
 
int GetBndCondIDToGlobalTraceID (const int i)
 Returns the global index of the boundary trace giving the index on the boundary expansion. More...
 
const Array< OneD, const int > & GetBndCondIDToGlobalTraceID ()
 
int GetNumGlobalDirBndCoeffs () const
 Returns the number of global Dirichlet boundary coefficients. More...
 
int GetNumLocalDirBndCoeffs () const
 Returns the number of local Dirichlet boundary coefficients. More...
 
int GetNumGlobalBndCoeffs () const
 Returns the total number of global boundary coefficients. More...
 
int GetNumLocalBndCoeffs () const
 Returns the total number of local boundary coefficients. More...
 
int GetNumLocalCoeffs () const
 Returns the total number of local coefficients. More...
 
int GetNumGlobalCoeffs () const
 Returns the total number of global coefficients. More...
 
bool GetSingularSystem () const
 Retrieves if the system is singular (true) or not (false) More...
 
void GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const
 
void GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
 
void GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc, int offset) const
 
void GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset, bool UseComm=true) const
 
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool UseComm=true) const
 
void LocalToLocalBnd (const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locbnd) const
 
void LocalToLocalInt (const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locint) const
 
void LocalBndToLocal (const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) const
 
void LocalIntToLocal (const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) const
 
void AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
 
void AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
void AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const
 
void AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal) const
 
void UniversalAssembleBnd (NekVector< NekDouble > &pGlobal) const
 
void UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal, int offset) const
 
int GetFullSystemBandWidth () const
 
int GetNumNonDirVertexModes () const
 
int GetNumNonDirEdgeModes () const
 
int GetNumNonDirFaceModes () const
 
int GetNumDirEdges () const
 
int GetNumDirFaces () const
 
int GetNumNonDirEdges () const
 
int GetNumNonDirFaces () const
 
void PrintStats (std::ostream &out, std::string variable, bool printHeader=true) const
 
const Array< OneD, const int > & GetExtraDirEdges ()
 
std::shared_ptr< AssemblyMap > LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
 
int GetBndSystemBandWidth () const
 Returns the bandwidth of the boundary system. More...
 
int GetStaticCondLevel () const
 Returns the level of static condensation for this map. More...
 
int GetNumPatches () const
 Returns the number of patches in this static condensation level. More...
 
const Array< OneD, const unsigned int > & GetNumLocalBndCoeffsPerPatch ()
 Returns the number of local boundary coefficients in each patch. More...
 
const Array< OneD, const unsigned int > & GetNumLocalIntCoeffsPerPatch ()
 Returns the number of local interior coefficients in each patch. More...
 
const AssemblyMapSharedPtr GetNextLevelLocalToGlobalMap () const
 Returns the local to global mapping for the next level in the multi-level static condensation. More...
 
void SetNextLevelLocalToGlobalMap (AssemblyMapSharedPtr pNextLevelLocalToGlobalMap)
 
const PatchMapSharedPtr & GetPatchMapFromPrevLevel (void) const
 Returns the patch map from the previous level of the multi-level static condensation. More...
 
bool AtLastLevel () const
 Returns true if this is the last level in the multi-level static condensation. More...
 
GlobalSysSolnType GetGlobalSysSolnType () const
 Returns the method of solving global systems. More...
 
std::string GetPreconType () const
 
bool IsAbsoluteTolerance () const
 
int GetSuccessiveRHS () const
 
std::string GetLinSysIterSolver () const
 
int GetLowestStaticCondLevel () const
 
void PatchLocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void PatchGlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
void PatchAssemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
int CreateGraph (const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
 
void SetUpUniversalC0ContMap (const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
 
void CalculateFullSystemBandWidth ()
 Calculate the bandwith of the full matrix system. More...
 
int v_GetLocalToGlobalMap (const int i) const override
 
int v_GetGlobalToUniversalMap (const int i) const override
 
int v_GetGlobalToUniversalMapUnique (const int i) const override
 
const Array< OneD, const int > & v_GetLocalToGlobalMap () override
 
const Array< OneD, const int > & v_GetGlobalToUniversalMap () override
 
const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique () override
 
NekDouble v_GetLocalToGlobalSign (const int i) const override
 
const Array< OneD, NekDouble > & v_GetLocalToGlobalSign () const override
 
void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const override
 
void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const override
 
void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const override
 
void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const override
 
void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const override
 
void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const override
 
void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const override
 
int v_GetFullSystemBandWidth () const override
 
int v_GetNumNonDirVertexModes () const override
 
int v_GetNumNonDirEdgeModes () const override
 
int v_GetNumNonDirFaceModes () const override
 
int v_GetNumDirEdges () const override
 
int v_GetNumDirFaces () const override
 
int v_GetNumNonDirEdges () const override
 
int v_GetNumNonDirFaces () const override
 
const Array< OneD, const int > & v_GetExtraDirEdges () override
 
AssemblyMapSharedPtr v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType) override
 Construct an AssemblyMapCG object which corresponds to the linear space of the current object. More...
 
- Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMap
void CalculateBndSystemBandWidth ()
 Calculates the bandwidth of the boundary system. More...
 
void GlobalToLocalBndWithoutSign (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc)
 
virtual int v_GetLocalToGlobalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMapUnique (const int i) const
 
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap ()
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap ()
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique ()
 
virtual NekDouble v_GetLocalToGlobalSign (const int i) const
 
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign () const
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
 
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
virtual void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
 
virtual void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
virtual void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
 
virtual int v_GetFullSystemBandWidth () const
 
virtual int v_GetNumNonDirVertexModes () const
 
virtual int v_GetNumNonDirEdgeModes () const
 
virtual int v_GetNumNonDirFaceModes () const
 
virtual int v_GetNumDirEdges () const
 
virtual int v_GetNumDirFaces () const
 
virtual int v_GetNumNonDirEdges () const
 
virtual int v_GetNumNonDirFaces () const
 
virtual const Array< OneD, const int > & v_GetExtraDirEdges ()
 
virtual std::shared_ptr< AssemblyMap > v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
 Generate a linear space mapping from existing mapping. More...
 
- Protected Attributes inherited from Nektar::MultiRegions::AssemblyMapCG
Array< OneD, int > m_localToGlobalMap
 Integer map of local coeffs to global space. More...
 
Array< OneD, NekDouble > m_localToGlobalSign
 Integer sign of local coeffs to global space. More...
 
int m_fullSystemBandWidth
 Bandwith of the full matrix system (no static condensation). More...
 
Array< OneD, int > m_globalToUniversalMap
 Integer map of process coeffs to universal space. More...
 
Array< OneD, int > m_globalToUniversalMapUnique
 Integer map of unique process coeffs to universal space (signed) More...
 
int m_numNonDirVertexModes
 Number of non Dirichlet vertex modes. More...
 
int m_numNonDirEdgeModes
 Number of non Dirichlet edge modes. More...
 
int m_numNonDirFaceModes
 Number of non Dirichlet face modes. More...
 
int m_numDirEdges
 Number of Dirichlet edges. More...
 
int m_numDirFaces
 Number of Dirichlet faces. More...
 
int m_numNonDirEdges
 Number of Dirichlet edges. More...
 
int m_numNonDirFaces
 Number of Dirichlet faces. More...
 
int m_numLocalBndCondCoeffs
 Number of local boundary condition coefficients. More...
 
Array< OneD, int > m_extraDirEdges
 Extra dirichlet edges in parallel. More...
 
int m_numLocDirBndCondDofs
 Number of local boundary condition degrees of freedom. More...
 
int m_maxStaticCondLevel
 Maximum static condensation level. More...
 
std::set< ExtraDirDof > m_copyLocalDirDofs
 Set indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. More...
 
std::set< int > m_parallelDirBndSign
 Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change. More...
 
- Protected Attributes inherited from Nektar::MultiRegions::AssemblyMap
LibUtilities::SessionReaderSharedPtr m_session
 Session object. More...
 
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
std::string m_variable
 Variable string identifier. More...
 
size_t m_hash
 Hash for map. More...
 
int m_numLocalBndCoeffs
 Number of local boundary coefficients. More...
 
int m_numGlobalBndCoeffs
 Total number of global boundary coefficients. More...
 
int m_numLocalDirBndCoeffs
 Number of Local Dirichlet Boundary Coefficients. More...
 
int m_numGlobalDirBndCoeffs
 Number of Global Dirichlet Boundary Coefficients. More...
 
bool m_systemSingular
 Flag indicating if the system is singular or not. More...
 
int m_numLocalCoeffs
 Total number of local coefficients. More...
 
int m_numGlobalCoeffs
 Total number of global coefficients. More...
 
bool m_signChange
 Flag indicating if modes require sign reversal. More...
 
Array< OneD, int > m_localToGlobalBndMap
 Integer map of local coeffs to global Boundary Dofs. More...
 
Array< OneD, NekDouble > m_localToGlobalBndSign
 Integer sign of local boundary coeffs to global space. More...
 
Array< OneD, int > m_localToLocalBndMap
 Integer map of local boundary coeffs to local boundary system numbering. More...
 
Array< OneD, int > m_localToLocalIntMap
 Integer map of local boundary coeffs to local interior system numbering. More...
 
Array< OneD, int > m_bndCondCoeffsToLocalCoeffsMap
 Integer map of bnd cond coeffs to local coefficients. More...
 
Array< OneD, NekDouble > m_bndCondCoeffsToLocalCoeffsSign
 Integer map of sign of bnd cond coeffs to local coefficients. More...
 
Array< OneD, int > m_bndCondCoeffsToLocalTraceMap
 Integer map of bnd cond coeff to local trace coeff. More...
 
Array< OneD, int > m_bndCondIDToGlobalTraceID
 Integer map of bnd cond trace number to global trace number. More...
 
Array< OneD, int > m_globalToUniversalBndMap
 Integer map of process coeffs to universal space. More...
 
Array< OneD, int > m_globalToUniversalBndMapUnique
 Integer map of unique process coeffs to universal space (signed) More...
 
GlobalSysSolnType m_solnType
 The solution type of the global system. More...
 
int m_bndSystemBandWidth
 The bandwith of the global bnd system. More...
 
std::string m_preconType
 Type type of preconditioner to use in iterative solver. More...
 
NekDouble m_iterativeTolerance
 Tolerance for iterative solver. More...
 
bool m_isAbsoluteTolerance
 
int m_successiveRHS
 sucessive RHS for iterative solver More...
 
std::string m_linSysIterSolver
 Iterative solver: Conjugate Gradient, GMRES. More...
 
Gs::gs_data * m_gsh
 
Gs::gs_data * m_bndGsh
 
Gs::gs_data * m_dirBndGsh
 gs gather communication to impose Dirhichlet BCs. More...
 
int m_staticCondLevel
 The level of recursion in the case of multi-level static condensation. More...
 
int m_numPatches
 The number of patches (~elements) in the current level. More...
 
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
 The number of bnd dofs per patch. More...
 
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
 The number of int dofs per patch. More...
 
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
 Map from the patches of the previous level to the patches of the current level. More...
 
int m_lowestStaticCondLevel
 Lowest static condensation level. More...
 

Detailed Description

Definition at line 45 of file CoupledLocalToGlobalC0ContMap.h.

Constructor & Destructor Documentation

◆ CoupledLocalToGlobalC0ContMap()

Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const SpatialDomains::MeshGraphSharedPtr &  graph,
const SpatialDomains::BoundaryConditionsSharedPtr &  boundaryConditions,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const MultiRegions::ExpListSharedPtr &  pressure,
const size_t  nz_loc,
const bool  CheckforSingularSys = true 
)

This is an vector extension of MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the Linearised Navier Stokes problem

STEP 1: Wrap boundary conditions vector in an array (since routine is set up for multiple fields) and call the graph re-odering subroutine to obtain the reordered values

@TODO: Fix this so that we can extract normals from edges

STEP 2a: Set the mean pressure modes to edges depending on type of direct solver technique;

STEP 2: Count out the number of Dirichlet vertices and edges first

STEP 3: Set up an array which contains the offset information of the different graph vertices.

This basically means to identify how many global degrees of freedom the individual graph vertices correspond. Obviously, the graph vertices corresponding to the mesh-vertices account for a single global DOF. However, the graph vertices corresponding to the element edges correspond to 2*(N-2) global DOF where N is equal to the number of boundary modes on this edge.

STEP 4: Now, all ingredients are ready to set up the actual local to global mapping.

The remainder of the map consists of the element-interior degrees of freedom. This leads to the block-diagonal submatrix as each element-interior mode is globally orthogonal to modes in all other elements.

STEP 5: The boundary condition mapping is generated from the same vertex renumbering and fill in a unique interior map.

Definition at line 51 of file CoupledLocalToGlobalC0ContMap.cpp.

59 : AssemblyMapCG(pSession, fields[0]->GetComm())
60{
61 size_t i, j, k, n;
62 size_t cnt = 0, offset = 0;
63 size_t meshVertId;
64 size_t meshEdgeId;
65 size_t bndEdgeCnt;
66 size_t globalId;
67 size_t nEdgeCoeffs;
68 size_t nEdgeInteriorCoeffs;
69 int firstNonDirGraphVertId;
70 size_t nLocBndCondDofs = 0;
71 size_t nLocDirBndCondDofs = 0;
72 int nExtraDirichlet = 0;
73 LocalRegions::Expansion2DSharedPtr locExpansion;
74 LocalRegions::SegExpSharedPtr bndSegExp;
75 LibUtilities::BasisType bType;
76 StdRegions::Orientation edgeOrient;
77 Array<OneD, unsigned int> edgeInteriorMap;
78 Array<OneD, int> edgeInteriorSign;
79 size_t nvel = fields.size();
80
81 const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
82 int id, diff;
83 size_t nel = fields[0]->GetNumElmts();
84
85 MultiRegions::PeriodicMap periodicVerts;
86 MultiRegions::PeriodicMap periodicEdges;
87 MultiRegions::PeriodicMap periodicFaces;
88 std::vector<std::map<int, int>> ReorderedGraphVertId(3);
89 MultiRegions::BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
90 int staticCondLevel = 0;
91
92 if (CheckforSingularSys) // all singularity checking by setting flag to true
93 {
94 m_systemSingular = true;
95 }
96 else // Turn off singular checking by setting flag to false
97 {
98 m_systemSingular = false;
99 }
100
101 /**
102 * STEP 1: Wrap boundary conditions vector in an array
103 * (since routine is set up for multiple fields) and call
104 * the graph re-odering subroutine to obtain the reordered
105 * values
106 */
107
108 // Obtain any periodic information and allocate default mapping array
109 fields[0]->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
110
111 const Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp =
112 fields[0]->GetBndCondExpansions();
113
114 Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr>>
115 bndConditionsVec(nvel);
116 for (i = 0; i < nvel; ++i)
117 {
118 bndConditionsVec[i] = fields[i]->GetBndConditions();
119 }
120
121 std::map<int, int> IsDirVertDof;
122 std::map<int, int> IsDirEdgeDof;
123
124 SpatialDomains::Geometry1DSharedPtr g;
125 for (j = 0; j < bndCondExp.size(); ++j)
126 {
127 std::map<int, int> BndExpVids;
128 // collect unique list of vertex ids for this expansion
129 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
130 {
131 g = bndCondExp[j]
132 ->GetExp(k)
133 ->as<LocalRegions::Expansion1D>()
134 ->GetGeom1D();
135 BndExpVids[g->GetVid(0)] = g->GetVid(0);
136 BndExpVids[g->GetVid(1)] = g->GetVid(1);
137 }
138
139 for (i = 0; i < nvel; ++i)
140 {
141 if (bndConditionsVec[i][j]->GetBoundaryConditionType() ==
142 SpatialDomains::eDirichlet)
143 {
144 // set number of Dirichlet conditions along edge
145 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
146 {
147 IsDirEdgeDof[bndCondExp[j]
148 ->GetExp(k)
149 ->as<LocalRegions::Expansion1D>()
150 ->GetGeom1D()
151 ->GetGlobalID()] += 1;
152 }
153
154 // Set number of Dirichlet conditions at vertices
155 // with a clamp on its maximum value being nvel to
156 // handle corners between expansions
157 for (auto &mapIt : BndExpVids)
158 {
159 id = IsDirVertDof[mapIt.second] + 1;
160 IsDirVertDof[mapIt.second] = (id > (int)nvel) ? nvel : id;
161 }
162 }
163 else
164 {
165 // Check to see that edge normals have non-zero
166 // component in this direction since otherwise
167 // also can be singular.
168 /// @TODO: Fix this so that we can extract normals from edges
169 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
170 {
171 Array<OneD, Array<OneD, NekDouble>> locnorm;
172 LocalRegions::Expansion1DSharedPtr loc_exp =
173 bndCondExp[j]
174 ->GetExp(k)
175 ->as<LocalRegions::Expansion1D>();
176 locnorm =
177 loc_exp->GetLeftAdjacentElementExp()->GetTraceNormal(
178 loc_exp->GetLeftAdjacentElementTrace());
179
180 size_t ndir = locnorm.size();
181 if (i < ndir) // account for Fourier version where n can be
182 // larger then ndir
183 {
184 for (size_t l = 0; l < locnorm[0].size(); ++l)
185 {
186 if (fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
187 {
188 m_systemSingular = false;
189 break;
190 }
191 }
192 }
193 if (m_systemSingular == false)
194 {
195 break;
196 }
197 }
198 }
199 }
200 }
201
202 Array<OneD, std::map<int, int>> Dofs(2);
203
204 Array<OneD, int> AddMeanPressureToEdgeId(nel, -1);
205 size_t edgeId, vertId;
206
207 // special case of singular problem - need to fix one pressure
208 // dof to a dirichlet edge. Since we attached pressure dof to
209 // last velocity component of edge need to make sure this
210 // component is Dirichlet
211 if (m_systemSingular)
212 {
213 id = -1;
214 for (i = 0; i < bndConditionsVec[0].size(); ++i)
215 {
216 if (bndConditionsVec[nvel - 1][i]->GetBoundaryConditionType() ==
217 SpatialDomains::eDirichlet)
218 {
219 id = bndCondExp[i]
220 ->GetExp(0)
221 ->as<LocalRegions::Expansion1D>()
222 ->GetGeom1D()
223 ->GetGlobalID();
224 break;
225 }
226 }
227
228 ASSERTL0(id != -1, " Did not find an edge to attach singular pressure "
229 "degree of freedom");
230
231 // determine element with this edge id. There may be a
232 // more direct way of getting element from spatialDomains
233 for (i = 0; i < nel; ++i)
234 {
235 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
236 {
237 edgeId = (locExpVector[i]
238 ->as<LocalRegions::Expansion2D>()
239 ->GetGeom2D())
240 ->GetEid(j);
241
242 if ((int)edgeId == id)
243 {
244 AddMeanPressureToEdgeId[i] = id;
245 break;
246 }
247 }
248
249 if (AddMeanPressureToEdgeId[i] != -1)
250 {
251 break;
252 }
253 }
254 }
255
256 for (i = 0; i < nel; ++i)
257 {
258 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
259 {
260 vertId =
261 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
262 ->GetVid(j);
263 if (Dofs[0].count(vertId) == 0)
264 {
265 Dofs[0][vertId] = nvel * nz_loc;
266
267 // Adjust for a Dirichlet boundary condition to give number to
268 // be solved
269 if (IsDirVertDof.count(vertId) != 0)
270 {
271 Dofs[0][vertId] -= IsDirVertDof[vertId] * nz_loc;
272 }
273 }
274
275 edgeId =
276 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
277 ->GetEid(j);
278 if (Dofs[1].count(edgeId) == 0)
279 {
280 Dofs[1][edgeId] =
281 nvel * (locExpVector[i]->GetTraceNcoeffs(j) - 2) * nz_loc;
282 }
283
284 // Adjust for Dirichlet boundary conditions to give number to be
285 // solved
286 if (IsDirEdgeDof.count(edgeId) != 0)
287 {
288 Dofs[1][edgeId] -= IsDirEdgeDof[edgeId] * nz_loc *
289 (locExpVector[i]->GetTraceNcoeffs(j) - 2);
290 }
291 }
292 }
293
294 std::set<int> extraDirVerts, extraDirEdges;
295
296 CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false, periodicVerts,
297 periodicEdges, periodicFaces, ReorderedGraphVertId,
298 bottomUpGraph, extraDirVerts, extraDirEdges,
299 firstNonDirGraphVertId, nExtraDirichlet, 4);
300 /*
301 SetUp2DGraphC0ContMap(*fields[0],
302 bndCondExp,
303 bndConditionsVec,
304 periodicVerts, periodicEdges,
305 Dofs, ReorderedGraphVertId,
306 firstNonDirGraphVertId, nExtraDirichlet,
307 bottomUpGraph, extraDir, false, 4);
308 */
309
310 /**
311 * STEP 2a: Set the mean pressure modes to edges depending on
312 * type of direct solver technique;
313 */
314
315 // determine which edge to add mean pressure dof based on
316 // ensuring that at least one pressure dof from an internal
317 // patch is associated with its boundary system
318 if (m_session->MatchSolverInfoAsEnum(
319 "GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
320 {
321
322 FindEdgeIdToAddMeanPressure(ReorderedGraphVertId, nel, locExpVector,
323 edgeId, vertId, firstNonDirGraphVertId,
324 IsDirEdgeDof, bottomUpGraph,
325 AddMeanPressureToEdgeId);
326 }
327
328 // Set unset elmts to non-Dirichlet edges.
329 // special case of singular problem - need to fix one
330 // pressure dof to a dirichlet edge
331 for (i = 0; i < nel; ++i)
332 {
333 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
334 {
335 edgeId =
336 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
337 ->GetEid(j);
338
339 if (IsDirEdgeDof.count(edgeId) == 0) // interior edge
340 {
341 // setup AddMeanPressureToEdgeId to decide where to
342 // put pressure
343 if (AddMeanPressureToEdgeId[i] == -1)
344 {
345 AddMeanPressureToEdgeId[i] = edgeId;
346 }
347 }
348 }
349 ASSERTL0((AddMeanPressureToEdgeId[i] != -1),
350 "Did not determine "
351 "an edge to attach mean pressure dof");
352 // Add the mean pressure degree of freedom to this edge
353 Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
354 }
355
356 std::map<int, int> pressureEdgeOffset;
357
358 /**
359 * STEP 2: Count out the number of Dirichlet vertices and edges first
360 */
361 for (i = 0; i < bndCondExp.size(); i++)
362 {
363 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
364 {
365 bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
366 for (k = 0; k < nvel; ++k)
367 {
368 if (bndConditionsVec[k][i]->GetBoundaryConditionType() ==
369 SpatialDomains::eDirichlet)
370 {
371 nLocDirBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
372 }
373
374 if (bndConditionsVec[k][i]->GetBoundaryConditionType() !=
375 SpatialDomains::ePeriodic)
376 {
377 nLocBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
378 }
379 }
380 }
381 }
382
383 if (m_systemSingular)
384 {
385 m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet + nz_loc;
386 }
387 else
388 {
389 m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
390 }
391
392 /**
393 * STEP 3: Set up an array which contains the offset information of
394 * the different graph vertices.
395 *
396 * This basically means to identify how many global degrees of
397 * freedom the individual graph vertices correspond. Obviously,
398 * the graph vertices corresponding to the mesh-vertices account
399 * for a single global DOF. However, the graph vertices
400 * corresponding to the element edges correspond to 2*(N-2) global DOF
401 * where N is equal to the number of boundary modes on this edge.
402 */
403 Array<OneD, int> graphVertOffset(
404 nvel * nz_loc *
405 (ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),
406 0);
407 graphVertOffset[0] = 0;
408
409 m_signChange = false;
410
411 for (i = 0; i < nel; ++i)
412 {
413 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
414
415 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
416 {
417 nEdgeCoeffs = locExpansion->GetTraceNcoeffs(j);
418 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
419 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
420
421 for (k = 0; k < nvel * nz_loc; ++k)
422 {
423 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
424 nz_loc +
425 k] = 1;
426 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] * nvel *
427 nz_loc +
428 k] = (nEdgeCoeffs - 2);
429 }
430
431 bType = locExpansion->GetBasisType(0);
432 // need a sign vector for modal expansions if nEdgeCoeffs >=4
433 if ((nEdgeCoeffs >= 4) && (bType == LibUtilities::eModified_A))
434 {
435 m_signChange = true;
436 }
437 }
438 }
439
440 // Add mean pressure modes;
441 for (i = 0; i < nel; ++i)
442 {
443 graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] +
444 1) *
445 nvel * nz_loc -
446 1] += nz_loc;
447 // graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc]
448 // += nz_loc;
449 }
450
451 // Negate the vertices and edges with only a partial
452 // Dirichlet conditon. Essentially we check to see if an edge
453 // has a mixed Dirichlet with Neumann/Robin Condition and if
454 // so negate the offset associated with this vertex.
455
456 std::map<int, int> DirVertChk;
457
458 for (i = 0; i < bndConditionsVec[0].size(); ++i)
459 {
460 cnt = 0;
461 for (j = 0; j < nvel; ++j)
462 {
463 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
464 SpatialDomains::eDirichlet)
465 {
466 cnt++;
467 }
468 }
469
470 // Case where partial Dirichlet boundary condition
471 if ((cnt > 0) && (cnt < nvel))
472 {
473 for (j = 0; j < nvel; ++j)
474 {
475 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
476 SpatialDomains::eDirichlet)
477 {
478 // negate graph offsets which should be
479 // Dirichlet conditions
480 for (k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
481 {
482 // vertices with mix condition;
483 id = bndCondExp[i]
484 ->GetExp(k)
485 ->as<LocalRegions::Expansion1D>()
486 ->GetGeom1D()
487 ->GetVid(0);
488 if (DirVertChk.count(id * nvel + j) == 0)
489 {
490 DirVertChk[id * nvel + j] = 1;
491 for (n = 0; n < nz_loc; ++n)
492 {
493 graphVertOffset[ReorderedGraphVertId[0][id] *
494 nvel * nz_loc +
495 j * nz_loc + n] *= -1;
496 }
497 }
498
499 id = bndCondExp[i]
500 ->GetExp(k)
501 ->as<LocalRegions::Expansion1D>()
502 ->GetGeom1D()
503 ->GetVid(1);
504 if (DirVertChk.count(id * nvel + j) == 0)
505 {
506 DirVertChk[id * nvel + j] = 1;
507 for (n = 0; n < nz_loc; ++n)
508 {
509 graphVertOffset[ReorderedGraphVertId[0][id] *
510 nvel * nz_loc +
511 j * nz_loc + n] *= -1;
512 }
513 }
514
515 // edges with mixed id;
516 id = bndCondExp[i]
517 ->GetExp(k)
518 ->as<LocalRegions::Expansion1D>()
519 ->GetGeom1D()
520 ->GetGlobalID();
521 for (n = 0; n < nz_loc; ++n)
522 {
523 graphVertOffset[ReorderedGraphVertId[1][id] * nvel *
524 nz_loc +
525 j * nz_loc + n] *= -1;
526 }
527 }
528 }
529 }
530 }
531 }
532
533 cnt = 0;
534 // assemble accumulative list of full Dirichlet values.
535 for (i = 0; i < firstNonDirGraphVertId * nvel * nz_loc; ++i)
536 {
537 diff = abs(graphVertOffset[i]);
538 graphVertOffset[i] = cnt;
539 cnt += diff;
540 }
541
542 // set Dirichlet values with negative values to Dirichlet value
543 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
544 ++i)
545 {
546 if (graphVertOffset[i] < 0)
547 {
548 diff = -graphVertOffset[i];
549 graphVertOffset[i] = -cnt;
550 cnt += diff;
551 }
552 }
553
554 // Accumulate all interior degrees of freedom with positive values
555 m_numGlobalDirBndCoeffs = cnt;
556
557 // offset values
558 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
559 ++i)
560 {
561 if (graphVertOffset[i] >= 0)
562 {
563 diff = graphVertOffset[i];
564 graphVertOffset[i] = cnt;
565 cnt += diff;
566 }
567 }
568
569 // Finally set negative entries (corresponding to Dirichlet
570 // values ) to be positive
571 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
572 ++i)
573 {
574 if (graphVertOffset[i] < 0)
575 {
576 graphVertOffset[i] = -graphVertOffset[i];
577 }
578 }
579
580 // Allocate the proper amount of space for the class-data and fill
581 // information that is already known
582 cnt = 0;
583 m_numLocalBndCoeffs = 0;
584 m_numLocalCoeffs = 0;
585
586 for (i = 0; i < nel; ++i)
587 {
588 m_numLocalBndCoeffs +=
589 nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
590 // add these coeffs up separately since
591 // pressure->GetNcoeffs can include the coefficient in
592 // multiple planes.
593 m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs() - 1) * nz_loc;
594 }
595
596 m_numLocalCoeffs += m_numLocalBndCoeffs;
597
598 m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs, -1);
599 m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs, -1);
600 m_bndCondIDToGlobalTraceID = Array<OneD, int>(nLocBndCondDofs, -1);
601
602 // Set default sign array.
603 m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs, 1.0);
604 m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs, 1.0);
605
606 m_staticCondLevel = staticCondLevel;
607 m_numPatches = nel;
608
609 m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
610 m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
611
612 for (i = 0; i < nel; ++i)
613 {
614 m_numLocalBndCoeffsPerPatch[i] =
615 (unsigned int)nz_loc *
616 (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
617 m_numLocalIntCoeffsPerPatch[i] =
618 (unsigned int)nz_loc * (pressure->GetExp(i)->GetNcoeffs() - 1);
619 }
620
621 // Set up local to local bnd and local int maps
622 m_localToLocalBndMap = Array<OneD, int>(m_numLocalBndCoeffs, -1);
623 m_localToLocalIntMap =
624 Array<OneD, int>(m_numLocalCoeffs - m_numLocalBndCoeffs, -1);
625
626 size_t bndcnt = 0;
627 size_t intcnt = 0;
628 cnt = 0;
629 for (i = 0; i < nel; ++i)
630 {
631 for (j = 0; j < nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs());
632 ++j)
633 {
634 m_localToLocalBndMap[bndcnt++] = cnt++;
635 }
636
637 for (n = 0; n < nz_loc; ++n)
638 {
639 m_localToLocalBndMap[bndcnt++] = cnt++;
640 for (j = 1; j < (size_t)pressure->GetExp(i)->GetNcoeffs(); ++j)
641 {
642 m_localToLocalIntMap[intcnt++] = cnt++;
643 }
644 }
645 }
646
647 /**
648 * STEP 4: Now, all ingredients are ready to set up the actual
649 * local to global mapping.
650 *
651 * The remainder of the map consists of the element-interior
652 * degrees of freedom. This leads to the block-diagonal submatrix
653 * as each element-interior mode is globally orthogonal to modes
654 * in all other elements.
655 */
656 cnt = 0;
657 size_t nv, velnbndry;
658 Array<OneD, unsigned int> bmap;
659
660 // Loop over all the elements in the domain in shuffled
661 // ordering (element type consistency)
662 for (i = 0; i < nel; ++i)
663 {
664 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
665
666 velnbndry = locExpansion->NumBndryCoeffs();
667
668 // Require an inverse ordering of the bmap system to store
669 // local numbering system. Therefore get hold of elemental
670 // bmap and set up an inverse map
671 std::map<int, int> inv_bmap;
672 locExpansion->GetBoundaryMap(bmap);
673 for (j = 0; j < bmap.size(); ++j)
674 {
675 inv_bmap[bmap[j]] = j;
676 }
677
678 // Loop over all edges (and vertices) of element i
679 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
680 {
681 nEdgeInteriorCoeffs = locExpansion->GetTraceNcoeffs(j) - 2;
682 edgeOrient = (locExpansion->GetGeom())->GetEorient(j);
683 meshEdgeId = (locExpansion->GetGeom())->GetEid(j);
684 meshVertId = (locExpansion->GetGeom())->GetVid(j);
685
686 auto pIt = periodicEdges.find(meshEdgeId);
687
688 // See if this edge is periodic. If it is, then we map all
689 // edges to the one with lowest ID, and align all
690 // coefficients to this edge orientation.
691 if (pIt != periodicEdges.end())
692 {
693 std::pair<int, StdRegions::Orientation> idOrient =
694 DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
695 pIt->second);
696 edgeOrient = idOrient.second;
697 }
698
699 locExpansion->GetTraceInteriorToElementMap(
700 j, edgeInteriorMap, edgeInteriorSign, edgeOrient);
701
702 // Set the global DOF for vertex j of element i
703 for (nv = 0; nv < nvel * nz_loc; ++nv)
704 {
705 m_localToGlobalMap[cnt + nv * velnbndry +
706 inv_bmap[locExpansion->GetVertexMap(j)]] =
707 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
708 nz_loc +
709 nv];
710
711 // Set the global DOF's for the interior modes of edge j
712 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
713 {
714 m_localToGlobalMap[cnt + nv * velnbndry +
715 inv_bmap[edgeInteriorMap[k]]] =
716 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] *
717 nvel * nz_loc +
718 nv] +
719 k;
720 }
721 }
722
723 // Fill the sign vector if required
724 if (m_signChange)
725 {
726 for (nv = 0; nv < nvel * nz_loc; ++nv)
727 {
728 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
729 {
730 m_localToGlobalSign[cnt + nv * velnbndry +
731 inv_bmap[edgeInteriorMap[k]]] =
732 (NekDouble)edgeInteriorSign[k];
733 }
734 }
735 }
736 }
737
738 // use difference between two edges of the
739 // AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
740 nEdgeInteriorCoeffs =
741 graphVertOffset[(ReorderedGraphVertId[1]
742 [AddMeanPressureToEdgeId[i]]) *
743 nvel * nz_loc +
744 1] -
745 graphVertOffset[(ReorderedGraphVertId[1]
746 [AddMeanPressureToEdgeId[i]]) *
747 nvel * nz_loc];
748
749 size_t psize = pressure->GetExp(i)->GetNcoeffs();
750 for (n = 0; n < nz_loc; ++n)
751 {
752 m_localToGlobalMap[cnt + nz_loc * nvel * velnbndry + n * psize] =
753 graphVertOffset
754 [(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] + 1) *
755 nvel * nz_loc -
756 1] +
757 nEdgeInteriorCoeffs +
758 pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
759
760 pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
761 }
762
763 cnt += (velnbndry * nvel + psize) * nz_loc;
764 }
765
766 // Set up the mapping for the boundary conditions
767 offset = cnt = 0;
768 for (nv = 0; nv < nvel; ++nv)
769 {
770 for (i = 0; i < bndCondExp.size(); i++)
771 {
772 if (bndConditionsVec[nv][i]->GetBoundaryConditionType() ==
773 SpatialDomains::ePeriodic)
774 {
775 continue;
776 }
777
778 for (n = 0; n < nz_loc; ++n)
779 {
780 int ncoeffcnt = 0;
781 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
782 {
783 bndSegExp =
784 bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
785
786 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
787 for (k = 0; k < 2; k++)
788 {
789 meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
790 m_bndCondIDToGlobalTraceID[cnt +
791 bndSegExp->GetVertexMap(k)] =
792 graphVertOffset[ReorderedGraphVertId[0]
793 [meshVertId] *
794 nvel * nz_loc +
795 nv * nz_loc + n];
796 }
797
798 meshEdgeId = (bndSegExp->GetGeom1D())->GetGlobalID();
799 bndEdgeCnt = 0;
800 nEdgeCoeffs = bndSegExp->GetNcoeffs();
801 for (k = 0; k < nEdgeCoeffs; k++)
802 {
803 if (m_bndCondIDToGlobalTraceID[cnt + k] == -1)
804 {
805 m_bndCondIDToGlobalTraceID[cnt + k] =
806 graphVertOffset
807 [ReorderedGraphVertId[1][meshEdgeId] *
808 nvel * nz_loc +
809 nv * nz_loc + n] +
810 bndEdgeCnt;
811 bndEdgeCnt++;
812 }
813 }
814 ncoeffcnt += nEdgeCoeffs;
815 }
816 // Note: Can not use bndCondExp[i]->GetNcoeffs()
817 // due to homogeneous extension not returning just
818 // the value per plane
819 offset += ncoeffcnt;
820 }
821 }
822 }
823
824 globalId = Vmath::Vmax(m_numLocalCoeffs, &m_localToGlobalMap[0], 1) + 1;
825 m_numGlobalBndCoeffs = globalId;
826
827 /**
828 * STEP 5: The boundary condition mapping is generated from the
829 * same vertex renumbering and fill in a unique interior map.
830 */
831 cnt = 0;
832 for (i = 0; i < (size_t)m_numLocalCoeffs; ++i)
833 {
834 if (m_localToGlobalMap[i] == -1)
835 {
836 m_localToGlobalMap[i] = globalId++;
837 }
838 else
839 {
840 if (m_signChange)
841 {
842 m_localToGlobalBndSign[cnt] = m_localToGlobalSign[i];
843 }
844 m_localToGlobalBndMap[cnt++] = m_localToGlobalMap[i];
845 }
846 }
847 m_numGlobalCoeffs = globalId;
848
849 // Set up the local to global map for the next level when using
850 // multi-level static condensation
851 if (m_session->MatchSolverInfoAsEnum(
852 "GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
853 {
854 if (m_staticCondLevel < (bottomUpGraph->GetNlevels() - 1))
855 {
856 Array<OneD, int> vwgts_perm(Dofs[0].size() + Dofs[1].size() -
857 firstNonDirGraphVertId);
858 for (i = 0; i < locExpVector.size(); ++i)
859 {
860 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
861 size_t nVert = locExpansion->GetNverts();
862 for (j = 0; j < nVert; ++j)
863 {
864 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
865 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
866
867 if (ReorderedGraphVertId[0][meshVertId] >=
868 firstNonDirGraphVertId)
869 {
870 vwgts_perm[ReorderedGraphVertId[0][meshVertId] -
871 firstNonDirGraphVertId] =
872 Dofs[0][meshVertId];
873 }
874
875 if (ReorderedGraphVertId[1][meshEdgeId] >=
876 firstNonDirGraphVertId)
877 {
878 vwgts_perm[ReorderedGraphVertId[1][meshEdgeId] -
879 firstNonDirGraphVertId] =
880 Dofs[1][meshEdgeId];
881 }
882 }
883 }
884
885 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
886
887 m_nextLevelLocalToGlobalMap =
888 MemoryManager<AssemblyMap>::AllocateSharedPtr(this,
889 bottomUpGraph);
890 }
891 }
892}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure
void FindEdgeIdToAddMeanPressure(std::vector< std::map< int, int > > &ReorderedGraphVertId, size_t &nel, const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId, size_t &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
Definition: CoupledLocalToGlobalC0ContMap.cpp:894
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
Nektar::MultiRegions::AssemblyMapCG::CreateGraph
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
Definition: AssemblyMapCG.cpp:77
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:106
Nektar::MultiRegions::AssemblyMapCG::AssemblyMapCG
AssemblyMapCG(const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &comm, const std::string variable="DefaultVar")
Default constructor.
Definition: AssemblyMapCG.cpp:69
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:108
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:361
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:375
Nektar::MultiRegions::AssemblyMap::m_localToLocalIntMap
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
Definition: AssemblyMap.h:384
Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:372
Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:380
Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffsPerPatch
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:427
Nektar::MultiRegions::AssemblyMap::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:330
Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:342
Nektar::MultiRegions::AssemblyMap::m_nextLevelLocalToGlobalMap
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:434
Nektar::MultiRegions::AssemblyMap::m_staticCondLevel
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:423
Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:346
Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:348
Nektar::MultiRegions::AssemblyMap::GetComm
LibUtilities::CommSharedPtr GetComm()
Retrieves the communicator.
Definition: AssemblyMap.cpp:693
Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:429
Nektar::MultiRegions::AssemblyMap::m_systemSingular
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:350
Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Definition: AssemblyMap.h:378
Nektar::MultiRegions::AssemblyMap::m_localToLocalBndMap
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
Definition: AssemblyMap.h:382
Nektar::MultiRegions::AssemblyMap::m_numPatches
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:425
Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:344
Nektar::MultiRegions::AssemblyMap::m_bndCondIDToGlobalTraceID
Array< OneD, int > m_bndCondIDToGlobalTraceID
Integer map of bnd cond trace number to global trace number.
Definition: AssemblyMap.h:392
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::SegExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:253
Nektar::LocalRegions::Expansion2DSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:46
Nektar::LocalRegions::Expansion1DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:50
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
Nektar::MultiRegions::BottomUpSubStructuredGraphSharedPtr
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
Definition: SubStructuredGraph.h:60
Nektar::MultiRegions::DeterminePeriodicEdgeOrientId
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
Definition: AssemblyMapCG.cpp:2232
Nektar::MultiRegions::PeriodicMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Definition: MultiRegions.hpp:175
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:49
Nektar::SpatialDomains::Geometry1DSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289

References tinysimd::abs(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::AssemblyMapCG::CreateGraph(), Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::LibUtilities::eModified_A, Nektar::SpatialDomains::ePeriodic, FindEdgeIdToAddMeanPressure(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Nektar::StdRegions::StdExpansion::GetVertexMap(), Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::AssemblyMap::m_bndCondIDToGlobalTraceID, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_localToLocalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToLocalIntMap, Nektar::MultiRegions::AssemblyMap::m_nextLevelLocalToGlobalMap, Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffsPerPatch, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch, Nektar::MultiRegions::AssemblyMap::m_numPatches, Nektar::MultiRegions::AssemblyMap::m_session, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::m_staticCondLevel, Nektar::MultiRegions::AssemblyMap::m_systemSingular, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), CG_Iterations::pressure, and Vmath::Vmax().

Member Function Documentation

◆ FindEdgeIdToAddMeanPressure()

void Nektar::CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure ( std::vector< std::map< int, int > > &  ReorderedGraphVertId,
size_t &  nel,
const LocalRegions::ExpansionVector &  locExpVector,
size_t &  edgeId,
size_t &  vertId,
int &  firstNonDirGraphVertId,
std::map< int, int > &  IsDirEdgeDof,
MultiRegions::BottomUpSubStructuredGraphSharedPtr &  bottomUpGraph,
Array< OneD, int > &  AddMeanPressureToEdgeId 
)

Definition at line 894 of file CoupledLocalToGlobalC0ContMap.cpp.

901{
902 size_t i, j, k;
903
904 // Make list of homogeneous graph edges to elmt mappings
905 Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(), 2, -1);
906 std::map<int, int> HomGraphEdgeIdToEdgeId;
907
908 for (i = 0; i < nel; ++i)
909 {
910 size_t nVert = locExpVector[i]->GetNverts();
911 for (j = 0; j < nVert; ++j)
912 {
913 edgeId =
914 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
915 ->GetEid(j);
916
917 // note second condition stops us using mixed boundary condition
918 if ((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId) &&
919 (IsDirEdgeDof.count(edgeId) == 0))
920 {
921 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId] -
922 firstNonDirGraphVertId] = edgeId;
923
924 if (EdgeIdToElmts[edgeId][0] == -1)
925 {
926 EdgeIdToElmts[edgeId][0] = i;
927 }
928 else
929 {
930 EdgeIdToElmts[edgeId][1] = i;
931 }
932 }
933 }
934 }
935
936 // Start at second to last level and find edge on boundary
937 // to attach element
938 size_t nlevels = bottomUpGraph->GetNlevels();
939
940 // determine a default edge to attach pressure modes to
941 // which is part of the inner solve;
942 int defedge = -1;
943
944 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
945 bottomUpGraph->GetInteriorBlocks(nlevels);
946 for (i = 0; i < bndgraphs.size(); ++i)
947 {
948 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
949 size_t nVert = bndgraphs[i]->GetNverts();
950
951 for (j = 0; j < nVert; ++j)
952 {
953 // find edge in graph vert list
954 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
955 {
956 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
957
958 if (defedge == -1)
959 {
960 defedge = edgeId;
961 break;
962 }
963 }
964 }
965 if (defedge != -1)
966 {
967 break;
968 }
969 }
970
971 for (size_t n = 1; n < nlevels; ++n)
972 {
973 // produce a map with a key that is the element id
974 // that contains which next level patch it belongs to
975 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
976 bottomUpGraph->GetInteriorBlocks(n + 1);
977
978 // Fill next level graph of adjacent elements and their level
979 std::map<int, int> ElmtInBndry;
980
981 for (i = 0; i < bndgraphs.size(); ++i)
982 {
983 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
984 size_t nVert = bndgraphs[i]->GetNverts();
985
986 for (j = 0; j < nVert; ++j)
987 {
988 // find edge in graph vert list
989 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
990 {
991 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
992
993 if (EdgeIdToElmts[edgeId][0] != -1)
994 {
995 ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
996 }
997 if (EdgeIdToElmts[edgeId][1] != -1)
998 {
999 ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
1000 }
1001 }
1002 }
1003 }
1004
1005 // Now search interior patches in this level for edges
1006 // that share the same element as a boundary edge and
1007 // assign this elmt that boundary edge
1008 std::vector<MultiRegions::SubGraphSharedPtr> intgraphs =
1009 bottomUpGraph->GetInteriorBlocks(n);
1010 for (i = 0; i < intgraphs.size(); ++i)
1011 {
1012 int GlobIdOffset = intgraphs[i]->GetIdOffset();
1013 bool SetEdge = false;
1014 int elmtid = 0;
1015 size_t nVert = intgraphs[i]->GetNverts();
1016 for (j = 0; j < nVert; ++j)
1017 {
1018 // Check to see if graph vert is an edge
1019 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1020 {
1021 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1022
1023 for (k = 0; k < 2; ++k)
1024 {
1025 // relevant edge id
1026 elmtid = EdgeIdToElmts[edgeId][k];
1027
1028 if (elmtid != -1)
1029 {
1030 auto mapIt = ElmtInBndry.find(elmtid);
1031
1032 if (mapIt != ElmtInBndry.end())
1033 {
1034 // now find a edge in the next level boundary
1035 // graph
1036 int GlobIdOffset1 =
1037 bndgraphs[mapIt->second]->GetIdOffset();
1038 for (int l = 0;
1039 l < bndgraphs[mapIt->second]->GetNverts();
1040 ++l)
1041 {
1042 // find edge in graph vert list
1043 if (HomGraphEdgeIdToEdgeId.count(
1044 GlobIdOffset1 + l) != 0)
1045 {
1046 AddMeanPressureToEdgeId[elmtid] =
1047 defedge;
1048
1049 SetEdge = true;
1050 break;
1051 }
1052 }
1053 }
1054 }
1055 }
1056 }
1057 }
1058
1059 // if we have failed to find matching edge in next
1060 // level patch boundary then set last found elmt
1061 // associated to this interior patch to the
1062 // default edget value
1063 if (SetEdge == false)
1064 {
1065 if (elmtid == -1) // find an elmtid in patch
1066 {
1067 size_t nVert = intgraphs[i]->GetNverts();
1068 for (j = 0; j < nVert; ++j)
1069 {
1070 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1071 {
1072 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1073 for (k = 0; k < 2; ++k)
1074 {
1075 // relevant edge id
1076 elmtid = EdgeIdToElmts[edgeId][k];
1077 if (elmtid != -1)
1078 {
1079 break;
1080 }
1081 }
1082 }
1083 if (elmtid != -1)
1084 {
1085 break;
1086 }
1087 }
1088 }
1089 if (AddMeanPressureToEdgeId[elmtid] == -1)
1090 {
1091 AddMeanPressureToEdgeId[elmtid] = defedge;
1092 }
1093 }
1094 }
1095 }
1096}

Referenced by CoupledLocalToGlobalC0ContMap().

Generated on Sun Jun 15 2025 07:52:36 for Nektar++ by doxygen 1.9.4