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 NekDoubleGetLocalToGlobalBndSign () 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< AssemblyMapLinearSpaceMap (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 PatchMapSharedPtrGetPatchMapFromPrevLevel (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
 
NekDouble GetIterativeTolerance () 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_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< 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 (NekVector< 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_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< 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 (NekVector< 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< AssemblyMapv_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, NekDoublem_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< ExtraDirDofm_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, NekDoublem_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, NekDoublem_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_datam_gsh
 
Gs::gs_datam_bndGsh
 
Gs::gs_datam_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 53 of file CoupledLocalToGlobalC0ContMap.cpp.

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

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 896 of file CoupledLocalToGlobalC0ContMap.cpp.

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

Referenced by CoupledLocalToGlobalC0ContMap().