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 (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 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::map< int, std::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, 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) 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...
 
std::map< int, std::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 53 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().

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

References Nektar::iterator.

Referenced by CoupledLocalToGlobalC0ContMap().

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