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...
 
virtual ~AssemblyMapCG ()
 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...
 
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...
 
PreconditionerType GetPreconType () const
 
NekDouble GetIterativeTolerance () const
 
int GetMaxIterations () 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...
 
virtual int v_GetLocalToGlobalMap (const int i) const override
 
virtual int v_GetGlobalToUniversalMap (const int i) const override
 
virtual int v_GetGlobalToUniversalMapUnique (const int i) const override
 
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap () override
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap () override
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique () override
 
virtual NekDouble v_GetLocalToGlobalSign (const int i) const override
 
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign () const override
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const override
 
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm) const override
 
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const override
 
virtual void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const override
 
virtual void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const override
 
virtual void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const override
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const override
 
virtual void v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const override
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const override
 
virtual int v_GetFullSystemBandWidth () const override
 
virtual int v_GetNumNonDirVertexModes () const override
 
virtual int v_GetNumNonDirEdgeModes () const override
 
virtual int v_GetNumNonDirFaceModes () const override
 
virtual int v_GetNumDirEdges () const override
 
virtual int v_GetNumDirFaces () const override
 
virtual int v_GetNumNonDirEdges () const override
 
virtual int v_GetNumNonDirFaces () const override
 
virtual const Array< OneD, const int > & v_GetExtraDirEdges () override
 
virtual 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)
 
- 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...
 
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...
 
PreconditionerType m_preconType
 Type type of preconditioner to use in iterative solver. More...
 
int m_maxIterations
 Maximum iterations for iterative solver. More...
 
NekDouble m_iterativeTolerance
 Tolerance for iterative solver. More...
 
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 53 of file CoupledLocalToGlobalC0ContMap.cpp.

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

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

Referenced by CoupledLocalToGlobalC0ContMap().

Generated on Sat Mar 25 2023 12:35:56 for Nektar++ by doxygen 1.9.1