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 int nz_loc, const bool CheeckForSingularSys=true)
 
void FindEdgeIdToAddMeanPressure (std::vector< std::map< int, int >> &ReorderedGraphVertId, int &nel, const LocalRegions::ExpansionVector &locExpVector, int &edgeId, int &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
 
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 AssemblyMapSharedPtr v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
 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 44 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 int  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 52 of file CoupledLocalToGlobalC0ContMap.cpp.

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

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::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,
int &  nel,
const LocalRegions::ExpansionVector &  locExpVector,
int &  edgeId,
int &  vertId,
int &  firstNonDirGraphVertId,
std::map< int, int > &  IsDirEdgeDof,
MultiRegions::BottomUpSubStructuredGraphSharedPtr &  bottomUpGraph,
Array< OneD, int > &  AddMeanPressureToEdgeId 
)

Definition at line 893 of file CoupledLocalToGlobalC0ContMap.cpp.

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

Referenced by CoupledLocalToGlobalC0ContMap().

Generated on Tue Aug 23 2022 07:27:18 for Nektar++ by doxygen 1.9.1