Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
Nektar::CoupledLocalToGlobalC0ContMap Class Reference

#include <CoupledLocalToGlobalC0ContMap.h>

Inheritance diagram for Nektar::CoupledLocalToGlobalC0ContMap:
Inheritance graph
[legend]
Collaboration diagram for Nektar::CoupledLocalToGlobalC0ContMap:
Collaboration graph
[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 (vector< map< int, int > > &ReorderedGraphVertId, int &nel, const LocalRegions::ExpansionVector &locExpVector, int &edgeId, int &vertId, int &firstNonDirGraphVertId, 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 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...
 
map< int, vector< ExtraDirDof > > & GetExtraDirDofs ()
 
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
 AssemblyMap ()
 Default constructor. More...
 
 AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, 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) const
 
void LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) 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
 
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...
 
int GetBndCondCoeffsToGlobalCoeffsMap (const int i)
 Retrieves the global index corresponding to a boundary expansion mode. More...
 
const Array< OneD, const int > & GetBndCondCoeffsToGlobalCoeffsMap ()
 Retrieves the global indices corresponding to the boundary expansion modes. More...
 
NekDouble GetBndCondCoeffsToGlobalCoeffsSign (const int i)
 Returns the modal sign associated with a given boundary expansion mode. More...
 
int GetBndCondTraceToGlobalTraceMap (const int i)
 Returns the global index of the boundary trace giving the index on the boundary expansion. More...
 
const Array< OneD, const int > & GetBndCondTraceToGlobalTraceMap ()
 
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 NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
 
void LocalBndToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const
 
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) 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) const
 
const Array< OneD, const int > & GetExtraDirEdges ()
 
boost::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
 
int GetLowestStaticCondLevel () 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, set< int > &extraDirVerts, 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) const
 
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) 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...
 
map< int, vector< ExtraDirDof > > m_extraDirDofs
 Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. 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 boundary coeffs to global space. More...
 
Array< OneD, NekDouble > m_localToGlobalBndSign
 Integer sign of local boundary coeffs to global space. More...
 
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
 Integer map of bnd cond coeffs to global coefficients. More...
 
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
 Integer map of bnd cond coeffs to global coefficients. More...
 
Array< OneD, int > m_bndCondTraceToGlobalTraceMap
 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...
 
Gs::gs_data * m_gsh
 
Gs::gs_data * m_bndGsh
 
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

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

: Fix this so that we can extract normals from edges

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

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

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

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

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

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

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

Definition at line 51 of file CoupledLocalToGlobalC0ContMap.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::AssemblyMapCG::CreateGraph(), Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, FindEdgeIdToAddMeanPressure(), Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::StdRegions::StdExpansion::GetVertexMap(), Nektar::iterator, Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign, 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(), and Vmath::Vmax().

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

Member Function Documentation

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

Definition at line 794 of file CoupledLocalToGlobalC0ContMap.cpp.

References Nektar::iterator.

Referenced by CoupledLocalToGlobalC0ContMap().

799 {
800 
801  int i,j,k;
802 
803  // Make list of homogeneous graph edges to elmt mappings
804  Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(),2,-1);
805  map<int,int> HomGraphEdgeIdToEdgeId;
806 
807  for(i = 0; i < nel; ++i)
808  {
809  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
810  {
811  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
812  ->GetGeom2D())->GetEid(j);
813 
814  // note second condition stops us using mixed boundary condition
815  if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
816  && (IsDirEdgeDof.count(edgeId) == 0))
817  {
818  HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
819 
820  if(EdgeIdToElmts[edgeId][0] == -1)
821  {
822  EdgeIdToElmts[edgeId][0] = i;
823  }
824  else
825  {
826  EdgeIdToElmts[edgeId][1] = i;
827  }
828  }
829  }
830  }
831 
832 
833  map<int,int>::iterator mapIt;
834 
835  // Start at second to last level and find edge on boundary
836  // to attach element
837  int nlevels = bottomUpGraph->GetNlevels();
838 
839  // determine a default edge to attach pressure modes to
840  // which is part of the inner solve;
841  int defedge = -1;
842 
843  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
844  for(i = 0; i < bndgraphs.size(); ++i)
845  {
846  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
847 
848  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
849  {
850  // find edge in graph vert list
851  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
852  {
853  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
854 
855  if(defedge == -1)
856  {
857  defedge = edgeId;
858  break;
859  }
860  }
861  }
862  if(defedge != -1)
863  {
864  break;
865  }
866  }
867 
868  for(int n = 1; n < nlevels; ++n)
869  {
870  // produce a map with a key that is the element id
871  // that contains which next level patch it belongs to
872  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
873 
874  // Fill next level graph of adjacent elements and their level
875  map<int,int> ElmtInBndry;
876 
877  for(i = 0; i < bndgraphs.size(); ++i)
878  {
879  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
880 
881  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
882  {
883  // find edge in graph vert list
884  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
885  {
886  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
887 
888  if(EdgeIdToElmts[edgeId][0] != -1)
889  {
890  ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
891  }
892  if(EdgeIdToElmts[edgeId][1] != -1)
893  {
894  ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
895  }
896  }
897  }
898  }
899 
900  // Now search interior patches in this level for edges
901  // that share the same element as a boundary edge and
902  // assign this elmt that boundary edge
903  vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
904  for(i = 0; i < intgraphs.size(); ++i)
905  {
906  int GlobIdOffset = intgraphs[i]->GetIdOffset();
907  bool SetEdge = false;
908  int elmtid = 0;
909  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
910  {
911  // Check to see if graph vert is an edge
912  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
913  {
914  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
915 
916  for(k = 0; k < 2; ++k)
917  {
918  // relevant edge id
919  elmtid = EdgeIdToElmts[edgeId][k];
920 
921  if(elmtid != -1)
922  {
923  mapIt = ElmtInBndry.find(elmtid);
924 
925  if(mapIt != ElmtInBndry.end())
926  {
927  // now find a edge in the next level boundary graph
928  int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
929  for(int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
930  {
931  // find edge in graph vert list
932  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
933  {
934  //June 2012: commenting this condition apparently
935  //solved the bug caused by the edge reordering procedure
936 
937  //if(AddMeanPressureToEdgeId[elmtid] == -1)
938  //{
939 
940  //AddMeanPressureToEdgeId[elmtid] = HomGraphEdgeIdToEdgeId[GlobIdOffset1+l];
941  AddMeanPressureToEdgeId[elmtid] = defedge;
942 
943  //}
944  SetEdge = true;
945  break;
946  }
947  }
948  }
949  }
950  }
951  }
952  }
953 
954 
955  // if we have failed to find matching edge in next
956  // level patch boundary then set last found elmt
957  // associated to this interior patch to the
958  // default edget value
959  if(SetEdge == false)
960  {
961  if(elmtid == -1) // find an elmtid in patch
962  {
963  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
964  {
965  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
966  {
967  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
968  for(k = 0; k < 2; ++k)
969  {
970  // relevant edge id
971  elmtid = EdgeIdToElmts[edgeId][k];
972  if(elmtid != -1)
973  {
974  break;
975  }
976  }
977  }
978  if(elmtid != -1)
979  {
980  break;
981  }
982  }
983  }
984  if(AddMeanPressureToEdgeId[elmtid] == -1)
985  {
986  AddMeanPressureToEdgeId[elmtid] = defedge;
987  }
988  }
989  }
990  }
991 
992 }
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Generated on Wed Mar 23 2016 20:14:28 for Nektar++ by   doxygen 1.8.8