Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | List of all members
Nektar::MultiRegions::AssemblyMapCG Class Reference

Constructs mappings for the C0 scalar continuous Galerkin formulation. More...

#include <AssemblyMapCG.h>

Inheritance diagram for Nektar::MultiRegions::AssemblyMapCG:
Inheritance graph
[legend]
Collaboration diagram for Nektar::MultiRegions::AssemblyMapCG:
Collaboration graph
[legend]

Public Member Functions

 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
 

Protected Member Functions

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

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...
 

Private Types

typedef Array< OneD, const
ExpListSharedPtr > 		BndCondExp
 
typedef Array< OneD, const
SpatialDomains::BoundaryConditionShPtr > 		BndCond
 

Detailed Description

Constructs mappings for the C0 scalar continuous Galerkin formulation.

Mappings are created for three possible global solution types:

These mappings are used by GlobalLinSys to generate the global system.

Definition at line 71 of file AssemblyMapCG.h.

Member Typedef Documentation

typedef Array<OneD, const SpatialDomains::BoundaryConditionShPtr> Nektar::MultiRegions::AssemblyMapCG::BndCond
private

Definition at line 75 of file AssemblyMapCG.h.

typedef Array<OneD, const ExpListSharedPtr> Nektar::MultiRegions::AssemblyMapCG::BndCondExp
private

Definition at line 73 of file AssemblyMapCG.h.

Constructor & Destructor Documentation

Nektar::MultiRegions::AssemblyMapCG::AssemblyMapCG ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const std::string  variable = "DefaultVar" 
)

Default constructor.

Definition at line 71 of file AssemblyMapCG.cpp.

References m_maxStaticCondLevel.

73  :
74  AssemblyMap(pSession,variable)
75  {
76  pSession->LoadParameter(
77  "MaxStaticCondLevel", m_maxStaticCondLevel, 100);
78  }
Nektar::MultiRegions::AssemblyMapCG::m_maxStaticCondLevel
int m_maxStaticCondLevel
Maximum static condensation level.
Definition: AssemblyMapCG.h:145
Nektar::MultiRegions::AssemblyMap::AssemblyMap
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:79
Nektar::MultiRegions::AssemblyMapCG::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.

STEP 6: 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.

Definition at line 1290 of file AssemblyMapCG.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::AssemblyMap::CalculateBndSystemBandWidth(), CalculateFullSystemBandWidth(), CreateGraph(), Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::MultiRegions::ePETScMultiLevelStaticCond, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Gs::Finalise(), Gs::Gather(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), Gs::gs_min, Gs::Init(), Nektar::iterator, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsMap, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsSign, Nektar::MultiRegions::AssemblyMap::m_comm, m_extraDirDofs, Nektar::MultiRegions::AssemblyMap::m_hash, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign, m_localToGlobalMap, 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, m_numLocalBndCondCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch, Nektar::MultiRegions::AssemblyMap::m_numPatches, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_staticCondLevel, Nektar::LibUtilities::ReduceSum, SetUpUniversalC0ContMap(), Nektar::MultiRegions::AssemblyMap::UniversalAssembleBnd(), and Vmath::Vmax().

1301  : AssemblyMap(pSession, variable)
1302  {
1303  int i, j, k, l;
1304  int cnt = 0;
1305  int intDofCnt;
1306  int meshVertId, meshEdgeId, meshFaceId;
1307  int globalId;
1308  int nEdgeInteriorCoeffs;
1309  int nFaceInteriorCoeffs;
1310  int firstNonDirGraphVertId;
1311  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1312  LocalRegions::ExpansionSharedPtr exp, bndExp;
1313  LibUtilities::BasisType bType;
1314  StdRegions::Orientation edgeOrient;
1315  StdRegions::Orientation faceOrient;
1316  Array<OneD, unsigned int> edgeInteriorMap;
1317  Array<OneD, int> edgeInteriorSign;
1318  Array<OneD, unsigned int> faceInteriorMap;
1319  Array<OneD, int> faceInteriorSign;
1320  PeriodicMap::const_iterator pIt;
1321 
1322  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1323 
1324  m_signChange = false;
1325 
1326  // Stores vertex, edge and face reordered vertices.
1327  DofGraph graph(3);
1328  DofGraph dofs(3);
1329 
1330  set<int> extraDirVerts, extraDirEdges;
1331  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
1332 
1333  // Construct list of number of degrees of freedom for each vertex,
1334  // edge and face.
1335  for (i = 0; i < locExpVector.size(); ++i)
1336  {
1337  exp = locExpVector[i];
1338 
1339  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
1340  {
1341  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1342  }
1343 
1344  for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
1345  {
1346  if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1347  {
1348  if (dofs[1][exp->GetGeom()->GetEid(j)] !=
1349  locExpVector[i]->GetEdgeNcoeffs(j)-2)
1350  {
1351  ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
1352  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B),
1353  "CG with variable order only available with modal expansion");
1354  }
1355  dofs[1][exp->GetGeom()->GetEid(j)] =
1356  min(dofs[1][exp->GetGeom()->GetEid(j)],
1357  locExpVector[i]->GetEdgeNcoeffs(j)-2);
1358  }
1359  else
1360  {
1361  dofs[1][exp->GetGeom()->GetEid(j)] =
1362  exp->GetEdgeNcoeffs(j) - 2;
1363  }
1364  }
1365 
1366  for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
1367  {
1368  if (dofs[2].count(exp->GetGeom()->GetFid(j)) > 0)
1369  {
1370  if (dofs[2][exp->GetGeom()->GetFid(j)] !=
1371  exp->GetFaceIntNcoeffs(j))
1372  {
1373  ASSERTL0( false,
1374  "CG with variable order not available in 3D");
1375  }
1376  dofs[2][exp->GetGeom()->GetFid(j)] =
1377  min(dofs[2][exp->GetGeom()->GetFid(j)],
1378  exp->GetFaceIntNcoeffs(j));
1379  }
1380  else
1381  {
1382  dofs[2][exp->GetGeom()->GetFid(j)] =
1383  exp->GetFaceIntNcoeffs(j);
1384  }
1385  }
1386  }
1387  // Now use information from all partitions to determine
1388  // the correct size
1389  map<int, int>::iterator dofIt;
1390  // edges
1391  Array<OneD, long> edgeId (dofs[1].size());
1392  Array<OneD, NekDouble> edgeDof (dofs[1].size());
1393  for(dofIt = dofs[1].begin(), i=0; dofIt != dofs[1].end(); dofIt++, i++)
1394  {
1395  edgeId[i] = dofIt->first;
1396  edgeDof[i] = (NekDouble) dofIt->second;
1397  }
1398  Gs::gs_data *tmp = Gs::Init(edgeId, vComm);
1399  Gs::Gather(edgeDof, Gs::gs_min, tmp);
1400  Gs::Finalise(tmp);
1401  for (i=0; i < dofs[1].size(); i++)
1402  {
1403  dofs[1][edgeId[i]] = (int) (edgeDof[i]+0.5);
1404  }
1405  // faces
1406  Array<OneD, long> faceId (dofs[2].size());
1407  Array<OneD, NekDouble> faceDof (dofs[2].size());
1408  for(dofIt = dofs[2].begin(), i=0; dofIt != dofs[2].end(); dofIt++, i++)
1409  {
1410  faceId[i] = dofIt->first;
1411  faceDof[i] = (NekDouble) dofIt->second;
1412  }
1413  Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
1414  Gs::Gather(faceDof, Gs::gs_min, tmp2);
1415  Gs::Finalise(tmp2);
1416  for (i=0; i < dofs[2].size(); i++)
1417  {
1418  dofs[2][faceId[i]] = (int) (faceDof[i]+0.5);
1419  }
1420 
1421  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1422 
1423  // Note that nExtraDirichlet is not used in the logic below; it just
1424  // needs to be set so that the coupled solver in
1425  // IncNavierStokesSolver can work.
1426  int nExtraDirichlet;
1427  int nGraphVerts =
1428  CreateGraph(locExp, bndCondExp, bndCondVec,
1429  checkIfSystemSingular, periodicVerts, periodicEdges,
1430  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1431  extraDirEdges, firstNonDirGraphVertId,
1432  nExtraDirichlet);
1433 
1434  /*
1435  * Set up an array which contains the offset information of the
1436  * different graph vertices.
1437  *
1438  * This basically means to identify to how many global degrees of
1439  * freedom the individual graph vertices correspond. Obviously,
1440  * the graph vertices corresponding to the mesh-vertices account
1441  * for a single global DOF. However, the graph vertices
1442  * corresponding to the element edges correspond to N-2 global DOF
1443  * where N is equal to the number of boundary modes on this edge.
1444  */
1445  Array<OneD, int> graphVertOffset(
1446  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1447 
1448  graphVertOffset[0] = 0;
1449 
1450  for(i = 0; i < locExpVector.size(); ++i)
1451  {
1452  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1453 
1454  for(j = 0; j < exp->GetNverts(); ++j)
1455  {
1456  meshVertId = exp->GetGeom()->GetVid(j);
1457  graphVertOffset[graph[0][meshVertId]+1] = 1;
1458  }
1459 
1460  for(j = 0; j < exp->GetNedges(); ++j)
1461  {
1462  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1463  meshEdgeId = exp->GetGeom()->GetEid(j);
1464  graphVertOffset[graph[1][meshEdgeId]+1]
1465  = dofs[1][meshEdgeId];
1466 
1467  bType = exp->GetEdgeBasisType(j);
1468 
1469  // need a sign vector for modal expansions if nEdgeCoeffs
1470  // >=3 (not 4 because of variable order case)
1471  if(nEdgeInteriorCoeffs &&
1472  (bType == LibUtilities::eModified_A ||
1473  bType == LibUtilities::eModified_B))
1474  {
1475  m_signChange = true;
1476  }
1477  }
1478 
1479  for(j = 0; j < exp->GetNfaces(); ++j)
1480  {
1481  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1482  meshFaceId = exp->GetGeom()->GetFid(j);
1483  graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
1484  }
1485  }
1486 
1487  for(i = 1; i < graphVertOffset.num_elements(); i++)
1488  {
1489  graphVertOffset[i] += graphVertOffset[i-1];
1490  }
1491 
1492  // Allocate the proper amount of space for the class-data
1493  m_numLocalCoeffs = numLocalCoeffs;
1494  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1495  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
1496  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
1497  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(m_numLocalBndCondCoeffs,-1);
1498  // If required, set up the sign-vector
1499  if(m_signChange)
1500  {
1501  m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
1502  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
1503  m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD,NekDouble>(m_numLocalBndCondCoeffs,1.0);
1504  }
1505 
1506  m_staticCondLevel = 0;
1507  m_numPatches = locExpVector.size();
1508  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1509  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1510  for(i = 0; i < m_numPatches; ++i)
1511  {
1512  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1513  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1514  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1515  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1516  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1517  }
1518 
1519  /**
1520  * STEP 6: Now, all ingredients are ready to set up the actual
1521  * local to global mapping.
1522  *
1523  * The remainder of the map consists of the element-interior
1524  * degrees of freedom. This leads to the block-diagonal submatrix
1525  * as each element-interior mode is globally orthogonal to modes
1526  * in all other elements.
1527  */
1528  cnt = 0;
1529 
1530  // Loop over all the elements in the domain
1531  for(i = 0; i < locExpVector.size(); ++i)
1532  {
1533  exp = locExpVector[i];
1534  cnt = locExp.GetCoeff_Offset(i);
1535  for(j = 0; j < exp->GetNverts(); ++j)
1536  {
1537  meshVertId = exp->GetGeom()->GetVid(j);
1538 
1539  // Set the global DOF for vertex j of element i
1540  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1541  graphVertOffset[graph[0][meshVertId]];
1542  }
1543 
1544  for(j = 0; j < exp->GetNedges(); ++j)
1545  {
1546  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1547  edgeOrient = exp->GetGeom()->GetEorient(j);
1548  meshEdgeId = exp->GetGeom()->GetEid(j);
1549 
1550  pIt = periodicEdges.find(meshEdgeId);
1551 
1552  // See if this edge is periodic. If it is, then we map all
1553  // edges to the one with lowest ID, and align all
1554  // coefficients to this edge orientation.
1555  if (pIt != periodicEdges.end())
1556  {
1557  pair<int, StdRegions::Orientation> idOrient =
1558  DeterminePeriodicEdgeOrientId(
1559  meshEdgeId, edgeOrient, pIt->second);
1560  edgeOrient = idOrient.second;
1561  }
1562 
1563  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1564 
1565  // Set the global DOF's for the interior modes of edge j
1566  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1567  {
1568  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1569  graphVertOffset[graph[1][meshEdgeId]]+k;
1570  }
1571  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1572  {
1573  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1574  graphVertOffset[graph[1][meshEdgeId]];
1575  }
1576 
1577  // Fill the sign vector if required
1578  if(m_signChange)
1579  {
1580  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1581  {
1582  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1583  }
1584  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1585  {
1586  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = 0.0;
1587  }
1588  }
1589  }
1590 
1591  for(j = 0; j < exp->GetNfaces(); ++j)
1592  {
1593  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1594  faceOrient = exp->GetGeom()->GetForient(j);
1595  meshFaceId = exp->GetGeom()->GetFid(j);
1596 
1597  pIt = periodicFaces.find(meshFaceId);
1598 
1599  if (pIt != periodicFaces.end() &&
1600  meshFaceId == min(meshFaceId, pIt->second[0].id))
1601  {
1602  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1603  }
1604 
1605  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1606 
1607  // Set the global DOF's for the interior modes of face j
1608  for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
1609  {
1610  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1611  graphVertOffset[graph[2][meshFaceId]]+k;
1612  }
1613  for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
1614  {
1615  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1616  graphVertOffset[graph[2][meshFaceId]];
1617  }
1618 
1619  if(m_signChange)
1620  {
1621  for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
1622  {
1623  m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
1624  }
1625  for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
1626  {
1627  m_localToGlobalSign[cnt+faceInteriorMap[k]] = 0;
1628  }
1629  }
1630 
1631  }
1632  }
1633 
1634  // Set up the mapping for the boundary conditions
1635  cnt = 0;
1636  int offset = 0;
1637  for(i = 0; i < bndCondExp.num_elements(); i++)
1638  {
1639  set<int> foundExtraVerts, foundExtraEdges;
1640  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1641  {
1642  bndExp = bndCondExp[i]->GetExp(j);
1643  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1644  for(k = 0; k < bndExp->GetNverts(); k++)
1645  {
1646  meshVertId = bndExp->GetGeom()->GetVid(k);
1647  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1648 
1649  if (bndConditions[i]->GetBoundaryConditionType() !=
1650  SpatialDomains::eDirichlet)
1651  {
1652  continue;
1653  }
1654 
1655  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1656  if (iter != extraDirVerts.end() &&
1657  foundExtraVerts.count(meshVertId) == 0)
1658  {
1659  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1660  bndExp->GetVertexMap(k);
1661  int gid = graphVertOffset[
1662  graph[0][meshVertId]];
1663  ExtraDirDof t(loc, gid, 1.0);
1664  m_extraDirDofs[i].push_back(t);
1665  foundExtraVerts.insert(meshVertId);
1666  }
1667  }
1668 
1669  for(k = 0; k < bndExp->GetNedges(); k++)
1670  {
1671  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1672  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1673  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1674 
1675  pIt = periodicEdges.find(meshEdgeId);
1676 
1677  // See if this edge is periodic. If it is, then we map
1678  // all edges to the one with lowest ID, and align all
1679  // coefficients to this edge orientation.
1680  if (pIt != periodicEdges.end())
1681  {
1682  pair<int, StdRegions::Orientation> idOrient =
1683  DeterminePeriodicEdgeOrientId(
1684  meshEdgeId, edgeOrient, pIt->second);
1685  edgeOrient = idOrient.second;
1686  }
1687 
1688  bndExp->GetEdgeInteriorMap(
1689  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1690 
1691  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1692  {
1693  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1694  graphVertOffset[graph[1][meshEdgeId]]+l;
1695  }
1696 
1697  // Fill the sign vector if required
1698  if(m_signChange)
1699  {
1700  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1701  {
1702  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1703  }
1704  }
1705 
1706  if (bndConditions[i]->GetBoundaryConditionType() !=
1707  SpatialDomains::eDirichlet)
1708  {
1709  continue;
1710  }
1711 
1712  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1713  if (iter != extraDirEdges.end() &&
1714  foundExtraEdges.count(meshEdgeId) == 0 &&
1715  nEdgeInteriorCoeffs > 0)
1716  {
1717  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1718  {
1719  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1720  edgeInteriorMap[l];
1721  int gid = graphVertOffset[
1722  graph[1][meshEdgeId]]+l;
1723  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1724  m_extraDirDofs[i].push_back(t);
1725  }
1726  foundExtraEdges.insert(meshEdgeId);
1727  }
1728  }
1729 
1730  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1731  intDofCnt = 0;
1732  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1733  {
1734  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1735  {
1736  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
1737  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1738  intDofCnt++;
1739  }
1740  }
1741  }
1742  offset += bndCondExp[i]->GetNcoeffs();
1743  }
1744 
1745  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1746  m_numGlobalBndCoeffs = globalId;
1747 
1748 
1749  /*
1750  * The boundary condition mapping is generated from the same vertex
1751  * renumbering.
1752  */
1753  cnt=0;
1754  for(i = 0; i < m_numLocalCoeffs; ++i)
1755  {
1756  if(m_localToGlobalMap[i] == -1)
1757  {
1758  m_localToGlobalMap[i] = globalId++;
1759  }
1760  else
1761  {
1762  if(m_signChange)
1763  {
1764  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
1765  }
1766  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
1767  }
1768  }
1769 
1770  m_numGlobalCoeffs = globalId;
1771 
1772  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1773 
1774  // Since we can now have multiple entries to m_extraDirDofs due to
1775  // periodic boundary conditions we make a call to work out the
1776  // multiplicity of all entries and invert (Need to be after
1777  // SetupUniversalC0ContMap)
1778  Array<OneD, NekDouble> valence(m_numGlobalBndCoeffs,0.0);
1779 
1780  // Fill in Dirichlet coefficients that are to be sent to other
1781  // processors with a value of 1
1782  map<int, vector<ExtraDirDof> >::iterator Tit;
1783 
1784  // Generate valence for extraDirDofs
1785  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1786  {
1787  for (i = 0; i < Tit->second.size(); ++i)
1788  {
1789  valence[Tit->second[i].get<1>()] = 1.0;
1790  }
1791  }
1792 
1793  UniversalAssembleBnd(valence);
1794 
1795  // Set third argument of tuple to inverse of valence.
1796  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1797  {
1798  for (i = 0; i < Tit->second.size(); ++i)
1799  {
1800  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1801  }
1802  }
1803 
1804  // Set up the local to global map for the next level when using
1805  // multi-level static condensation
1806  if ((m_solnType == eDirectMultiLevelStaticCond ||
1807  m_solnType == eIterativeMultiLevelStaticCond ||
1808  m_solnType == eXxtMultiLevelStaticCond ||
1809  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1810  {
1811  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1812  {
1813  Array<OneD, int> vwgts_perm(
1814  dofs[0].size() + dofs[1].size() + dofs[2].size()
1815  - firstNonDirGraphVertId);
1816 
1817  for (i = 0; i < locExpVector.size(); ++i)
1818  {
1819  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1820 
1821  for (j = 0; j < exp->GetNverts(); ++j)
1822  {
1823  meshVertId = exp->GetGeom()->GetVid(j);
1824 
1825  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1826  {
1827  vwgts_perm[graph[0][meshVertId] -
1828  firstNonDirGraphVertId] =
1829  dofs[0][meshVertId];
1830  }
1831  }
1832 
1833  for (j = 0; j < exp->GetNedges(); ++j)
1834  {
1835  meshEdgeId = exp->GetGeom()->GetEid(j);
1836 
1837  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
1838  {
1839  vwgts_perm[graph[1][meshEdgeId] -
1840  firstNonDirGraphVertId] =
1841  dofs[1][meshEdgeId];
1842  }
1843  }
1844 
1845  for (j = 0; j < exp->GetNfaces(); ++j)
1846  {
1847  meshFaceId = exp->GetGeom()->GetFid(j);
1848 
1849  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
1850  {
1851  vwgts_perm[graph[2][meshFaceId] -
1852  firstNonDirGraphVertId] =
1853  dofs[2][meshFaceId];
1854  }
1855  }
1856  }
1857 
1858  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1859  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
1860  AllocateSharedPtr(this, bottomUpGraph);
1861  }
1862  }
1863 
1864  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
1865  m_localToGlobalMap.end());
1866 
1867  // Add up hash values if parallel
1868  int hash = m_hash;
1869  vComm->AllReduce(hash, LibUtilities::ReduceSum);
1870  m_hash = hash;
1871 
1872  CalculateBndSystemBandWidth();
1873  CalculateFullSystemBandWidth();
1874  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
Nektar::MultiRegions::AssemblyMap::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
Gs::Gather
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:224
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Gs::Finalise
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:206
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::SetUpUniversalC0ContMap
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
Definition: AssemblyMapCG.cpp:1995
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
Gs::Init
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
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::MultiRegions::AssemblyMapCG::m_numLocalBndCondCoeffs
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
Definition: AssemblyMapCG.h:139
Nektar::MultiRegions::AssemblyMapCG::m_extraDirDofs
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
Definition: AssemblyMapCG.h:148
Nektar::MultiRegions::AssemblyMapCG::CalculateFullSystemBandWidth
void CalculateFullSystemBandWidth()
Calculate the bandwith of the full matrix system.
Definition: AssemblyMapCG.cpp:2340
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::MultiRegions::AssemblyMap::m_hash
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:309
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
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_solnType
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:363
Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsSign
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:354
Nektar::MultiRegions::ExtraDirDof
boost::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:54
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::MultiRegions::DeterminePeriodicFaceOrient
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:1952
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LocalRegions::ExpansionSharedPtr
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
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::CalculateBndSystemBandWidth
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Definition: AssemblyMap.cpp:456
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_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::DofGraph
std::vector< std::map< int, int > > DofGraph
Definition: AssemblyMapCG.h:56
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::MultiRegions::AssemblyMap::AssemblyMap
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:79
Gs::gs_data
Definition: GsLib.hpp:117
Nektar::MultiRegions::AssemblyMap::UniversalAssembleBnd
void UniversalAssembleBnd(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:1128
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:1904
Nektar::MultiRegions::AssemblyMap::m_numPatches
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:386
Nektar::MultiRegions::AssemblyMapCG::~AssemblyMapCG ( )
virtual

Destructor.

Definition at line 1879 of file AssemblyMapCG.cpp.

References Gs::Finalise(), Nektar::MultiRegions::AssemblyMap::m_bndGsh, and Nektar::MultiRegions::AssemblyMap::m_gsh.

1880  {
1881  Gs::Finalise(m_gsh);
1882  Gs::Finalise(m_bndGsh);
1883  }
Gs::Finalise
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:206

Member Function Documentation

void Nektar::MultiRegions::AssemblyMapCG::CalculateFullSystemBandWidth ( )
protected

Calculate the bandwith of the full matrix system.

The bandwidth calculated here corresponds to what is referred to as half-bandwidth. If the elements of the matrix are designated as a_ij, it corresponds to the maximum value of |i-j| for non-zero a_ij. As a result, the value also corresponds to the number of sub- or super-diagonals.

The bandwith can be calculated elementally as it corresponds to the maximal elemental bandwith (i.e. the maximal difference in global DOF index for every element).

We caluclate here the bandwith of the full global system.

Definition at line 2340 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth, m_localToGlobalMap, Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffsPerPatch, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch, and Nektar::MultiRegions::AssemblyMap::m_numPatches.

Referenced by AssemblyMapCG().

2341  {
2342  int i,j;
2343  int cnt = 0;
2344  int locSize;
2345  int maxId;
2346  int minId;
2347  int bwidth = -1;
2348  for(i = 0; i < m_numPatches; ++i)
2349  {
2350  locSize = m_numLocalBndCoeffsPerPatch[i]+m_numLocalIntCoeffsPerPatch[i];
2351  maxId = -1;
2352  minId = m_numLocalCoeffs+1;
2353  for(j = 0; j < locSize; j++)
2354  {
2355  if(m_localToGlobalMap[cnt+j] >= m_numGlobalDirBndCoeffs)
2356  {
2357  if(m_localToGlobalMap[cnt+j] > maxId)
2358  {
2359  maxId = m_localToGlobalMap[cnt+j];
2360  }
2361 
2362  if(m_localToGlobalMap[cnt+j] < minId)
2363  {
2364  minId = m_localToGlobalMap[cnt+j];
2365  }
2366  }
2367  }
2368  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2369 
2370  cnt+=locSize;
2371  }
2372 
2373  m_fullSystemBandWidth = bwidth;
2374  }
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
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_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::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:390
Nektar::MultiRegions::AssemblyMapCG::m_fullSystemBandWidth
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
Definition: AssemblyMapCG.h:119
Nektar::MultiRegions::AssemblyMap::m_numPatches
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:386
int Nektar::MultiRegions::AssemblyMapCG::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 
)
protected
  • Count verts, edges, face and add up edges and face sizes
  • Periodic vertices
  • Periodic edges
  • Periodic faces

STEP 4: Fill the #graph[0] and #graph[1] with the optimal ordering from boost.

Definition at line 80 of file AssemblyMapCG.cpp.

References ASSERTL0, ASSERTL1, Nektar::MultiRegions::CuthillMckeeReordering(), Nektar::MultiRegions::eDirectFullMatrix, Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::MultiRegions::eDirectStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::MultiRegions::eIterativeFull, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::MultiRegions::eIterativeStaticCond, Nektar::SpatialDomains::eNeumann, Nektar::MultiRegions::ePETScFullMatrix, Nektar::MultiRegions::ePETScMultiLevelStaticCond, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::eXxtFullMatrix, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Nektar::MultiRegions::eXxtStaticCond, Gs::Finalise(), Gs::Gather(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), Nektar::MultiRegions::GlobalSysSolnTypeMap, Gs::gs_add, Vmath::Imax(), Gs::Init(), Nektar::iterator, Nektar::MultiRegions::AssemblyMap::m_comm, m_extraDirEdges, Nektar::MultiRegions::AssemblyMap::m_lowestStaticCondLevel, m_numDirEdges, m_numDirFaces, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs, m_numLocalBndCondCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs, m_numNonDirEdgeModes, m_numNonDirEdges, m_numNonDirFaceModes, m_numNonDirFaces, m_numNonDirVertexModes, Nektar::MultiRegions::AssemblyMap::m_session, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_systemSingular, Nektar::MultiRegions::MultiLevelBisectionReordering(), Nektar::MultiRegions::NoReordering(), Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceMin, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

Referenced by AssemblyMapCG(), and Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap().

95  {
96  int graphVertId = 0;
97  int vMaxVertId = -1;
98  int i, j, k, l, cnt;
99  int meshVertId, meshEdgeId, meshFaceId;
100  int meshVertId2, meshEdgeId2;
101 
102  LocalRegions::ExpansionSharedPtr exp, bndExp;
103  const LocalRegions::ExpansionVector &locExpVector =
104  *(locExp.GetExp());
105  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
106  PeriodicMap::const_iterator pIt;
107 
108  m_numLocalBndCondCoeffs = 0;
109  m_systemSingular = checkIfSystemSingular;
110 
111  for(i = 0; i < bndCondExp.num_elements(); i++)
112  {
113  // Check to see if any value on boundary has Dirichlet value.
114  cnt = 0;
115  for(k = 0; k < bndConditions.num_elements(); ++k)
116  {
117  if (bndConditions[k][i]->GetBoundaryConditionType() ==
118  SpatialDomains::eDirichlet)
119  {
120  cnt++;
121  }
122  if (bndConditions[k][i]->GetBoundaryConditionType() !=
123  SpatialDomains::eNeumann)
124  {
125  m_systemSingular = false;
126  }
127  }
128 
129  // Find the maximum boundary vertex ID on this process. This is
130  // used later to pin a vertex if the system is singular.
131  for (j = 0; j < bndCondExp[i]->GetNumElmts(); ++j)
132  {
133  bndExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::Expansion>();
134  for (k = 0; k < bndExp->GetNverts(); ++k)
135  {
136  if (vMaxVertId < bndExp->GetGeom()->GetVid(k))
137  {
138  vMaxVertId = bndExp->GetGeom()->GetVid(k);
139  }
140  }
141 
142  }
143 
144  // If all boundaries are Dirichlet take out of mask
145  if(cnt == bndConditions.num_elements())
146  {
147  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
148  {
149  bndExp = bndCondExp[i]->GetExp(j);
150 
151  for (k = 0; k < bndExp->GetNverts(); k++)
152  {
153  meshVertId = bndExp->GetGeom()->GetVid(k);
154  if (graph[0].count(meshVertId) == 0)
155  {
156  graph[0][meshVertId] = graphVertId++;
157  }
158  }
159 
160  for (k = 0; k < bndExp->GetNedges(); k++)
161  {
162  meshEdgeId = bndExp->GetGeom()->GetEid(k);
163  if (graph[1].count(meshEdgeId) == 0)
164  {
165  graph[1][meshEdgeId] = graphVertId++;
166  }
167  }
168 
169  // Possibly not a face.
170  meshFaceId = bndExp->GetGeom()->GetGlobalID();
171  const int bndDim = bndExp->GetNumBases();
172  if (graph[bndDim].count(meshFaceId) == 0)
173  {
174  graph[bndDim][meshFaceId] = graphVertId++;
175  }
176  m_numLocalDirBndCoeffs += bndExp->GetNcoeffs();
177  }
178  }
179 
180  m_numLocalBndCondCoeffs += bndCondExp[i]->GetNcoeffs();
181  }
182 
183  // Number of dirichlet edges and faces (not considering periodic
184  // BCs)
185  m_numDirEdges = graph[1].size();
186  m_numDirFaces = graph[2].size();
187 
188  /*
189  * The purpose of this routine is to deal with those degrees of
190  * freedom that are Dirichlet, but do not have a local Dirichlet
191  * boundary condition expansion set.
192  *
193  * For example, in 2D, consider a triangulation of a square into two
194  * triangles. Now imagine one edge of the square is Dirichlet and
195  * the problem is run on two processors. On one processor, one
196  * triangle vertex is Dirichlet, but doesn't know this since the
197  * Dirichlet composite lives on the other processor.
198  *
199  * When the global linear system is solved therefore, there is an
200  * inconsistency that at best leads to an inaccurate answer or a
201  * divergence of the system.
202  *
203  * This routine identifies such cases for 2D, and also for 3D where
204  * e.g. edges may have the same problem (consider an extrusion of
205  * the case above, for example).
206  */
207 
208  // Collate information on Dirichlet vertices from all processes
209  int n = vComm->GetSize();
210  int p = vComm->GetRank();
211 
212  // At this point, graph only contains information from Dirichlet
213  // boundaries. Therefore make a global list of the vert and edge
214  // information on all processors.
215  Array<OneD, int> vertcounts (n, 0);
216  Array<OneD, int> vertoffsets(n, 0);
217  Array<OneD, int> edgecounts (n, 0);
218  Array<OneD, int> edgeoffsets(n, 0);
219  vertcounts[p] = graph[0].size();
220  edgecounts[p] = graph[1].size();
221  vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
222  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
223 
224  for (i = 1; i < n; ++i)
225  {
226  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
227  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
228  }
229 
230  int nTotVerts = Vmath::Vsum(n,vertcounts,1);
231  int nTotEdges = Vmath::Vsum(n,edgecounts,1);
232 
233  Array<OneD, int> vertlist(nTotVerts, 0);
234  Array<OneD, int> edgelist(nTotEdges, 0);
235  std::map<int, int>::iterator it;
236 
237  // construct list of global ids of global vertices
238  for (it = graph[0].begin(), i = 0;
239  it != graph[0].end();
240  ++it, ++i)
241  {
242  vertlist[vertoffsets[p] + i] = it->first;
243  }
244 
245  // construct list of global ids of global edges
246  for (it = graph[1].begin(), i = 0;
247  it != graph[1].end();
248  ++it, ++i)
249  {
250  edgelist[edgeoffsets[p] + i] = it->first;
251  }
252  vComm->AllReduce(vertlist, LibUtilities::ReduceSum);
253  vComm->AllReduce(edgelist, LibUtilities::ReduceSum);
254 
255  // Now we have a list of all Dirichlet vertices and edges on all
256  // processors.
257  nExtraDirichlet = 0;
258  map<int, int> extraDirVertIds, extraDirEdgeIds;
259 
260  // Ensure Dirchlet vertices are consistently recorded between
261  // processes (e.g. Dirichlet region meets Neumann region across a
262  // partition boundary requires vertex on partition to be Dirichlet).
263  //
264  // To do this we look over all elements and vertices in local
265  // partition and see if they match the values stored in the vertlist
266  // from othe processors and if so record the meshVertId/meshEdgeId
267  // and the processor it comes from.
268  for (i = 0; i < n; ++i)
269  {
270  if (i == p)
271  {
272  continue;
273  }
274 
275  for(j = 0; j < locExpVector.size(); j++)
276  {
277  exp = locExpVector[locExp.GetOffset_Elmt_Id(j)];
278 
279  for(k = 0; k < exp->GetNverts(); k++)
280  {
281  meshVertId = exp->GetGeom()->GetVid(k);
282  if(graph[0].count(meshVertId) == 0)
283  {
284  for (l = 0; l < vertcounts[i]; ++l)
285  {
286  if (vertlist[vertoffsets[i]+l] == meshVertId)
287  {
288  extraDirVertIds[meshVertId] = i;
289  graph[0][meshVertId] = graphVertId++;
290  nExtraDirichlet++;
291  }
292  }
293  }
294  }
295 
296  for(k = 0; k < exp->GetNedges(); k++)
297  {
298  meshEdgeId = exp->GetGeom()->GetEid(k);
299  if(graph[1].count(meshEdgeId) == 0)
300  {
301  for (l = 0; l < edgecounts[i]; ++l)
302  {
303  if (edgelist[edgeoffsets[i]+l] == meshEdgeId)
304  {
305  extraDirEdgeIds[meshEdgeId] = i;
306  graph[1][meshEdgeId] = graphVertId++;
307  nExtraDirichlet += exp->GetEdgeNcoeffs(k)-2;
308  }
309  }
310  }
311  }
312  }
313  }
314 
315  // Low Energy preconditioner needs to know how many extra Dirichlet
316  // edges are on this process so store map in array.
317  map<int, int>::const_iterator mapConstIt;
318  m_extraDirEdges = Array<OneD, int>(extraDirEdgeIds.size(), -1);
319  for (mapConstIt = extraDirEdgeIds.begin(), i = 0;
320  mapConstIt != extraDirEdgeIds.end(); mapConstIt++)
321  {
322  meshEdgeId = mapConstIt->first;
323  m_extraDirEdges[i++] = meshEdgeId;
324  }
325 
326  // Now we have a list of all vertices and edges that are Dirichlet
327  // and not defined on the local partition as well as which processor
328  // they are stored on.
329  //
330  // Make a full list of all such entities on all processors and which
331  // processor they belong to.
332  for (i = 0; i < n; ++i)
333  {
334  vertcounts [i] = 0;
335  vertoffsets[i] = 0;
336  edgecounts [i] = 0;
337  edgeoffsets[i] = 0;
338  }
339 
340  vertcounts[p] = extraDirVertIds.size();
341  edgecounts[p] = extraDirEdgeIds.size();
342  vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
343  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
344  nTotVerts = Vmath::Vsum(n, vertcounts, 1);
345  nTotEdges = Vmath::Vsum(n, edgecounts, 1);
346 
347  vertoffsets[0] = edgeoffsets[0] = 0;
348 
349  for (i = 1; i < n; ++i)
350  {
351  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
352  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
353  }
354 
355  Array<OneD, int> vertids (nTotVerts, 0);
356  Array<OneD, int> edgeids (nTotEdges, 0);
357  Array<OneD, int> vertprocs(nTotVerts, 0);
358  Array<OneD, int> edgeprocs(nTotEdges, 0);
359 
360  for (it = extraDirVertIds.begin(), i = 0;
361  it != extraDirVertIds.end(); ++it, ++i)
362  {
363  vertids [vertoffsets[p]+i] = it->first;
364  vertprocs[vertoffsets[p]+i] = it->second;
365  }
366 
367  for (it = extraDirEdgeIds.begin(), i = 0;
368  it != extraDirEdgeIds.end(); ++it, ++i)
369  {
370  edgeids [edgeoffsets[p]+i] = it->first;
371  edgeprocs[edgeoffsets[p]+i] = it->second;
372  }
373 
374  vComm->AllReduce(vertids, LibUtilities::ReduceSum);
375  vComm->AllReduce(vertprocs, LibUtilities::ReduceSum);
376  vComm->AllReduce(edgeids, LibUtilities::ReduceSum);
377  vComm->AllReduce(edgeprocs, LibUtilities::ReduceSum);
378 
379  // Set up list of vertices that need to be shared to other
380  // partitions
381  for (i = 0; i < nTotVerts; ++i)
382  {
383  if (vComm->GetRank() == vertprocs[i])
384  {
385  extraDirVerts.insert(vertids[i]);
386  }
387  }
388 
389  // Set up list of edges that need to be shared to other partitions
390  for (i = 0; i < nTotEdges; ++i)
391  {
392  if (vComm->GetRank() == edgeprocs[i])
393  {
394  extraDirEdges.insert(edgeids[i]);
395  }
396  }
397 
398  // Check between processes if the whole system is singular
399  int s = m_systemSingular ? 1 : 0;
400  vComm->AllReduce(s, LibUtilities::ReduceMin);
401  m_systemSingular = s == 1 ? true : false;
402 
403  // Find the minimum boundary vertex ID on each process
404  Array<OneD, int> bcminvertid(n, 0);
405  bcminvertid[p] = vMaxVertId;
406  vComm->AllReduce(bcminvertid, LibUtilities::ReduceMax);
407 
408  // Find the process rank with the minimum boundary vertex ID
409  int maxIdx = Vmath::Imax(n, bcminvertid, 1);
410 
411  // If the system is singular, the process with the maximum number of
412  // BCs will set a Dirichlet vertex to make system non-singular.
413  // Note: we find the process with maximum boundary regions to ensure
414  // we do not try to set a Dirichlet vertex on a partition with no
415  // intersection with the boundary.
416  meshVertId = 0;
417 
418  if (m_systemSingular && checkIfSystemSingular && maxIdx == p)
419  {
420  if (m_session->DefinesParameter("SingularVertex"))
421  {
422  m_session->LoadParameter("SingularVertex", meshVertId);
423  }
424  else if (vMaxVertId == -1)
425  {
426  // All boundaries are periodic.
427  meshVertId = locExpVector[0]->GetGeom()->GetVid(0);
428  }
429  else
430  {
431  // Set pinned vertex to that with minimum vertex ID to
432  // ensure consistency in parallel.
433  meshVertId = bcminvertid[p];
434  }
435 
436  if (graph[0].count(meshVertId) == 0)
437  {
438  graph[0][meshVertId] = graphVertId++;
439  }
440  }
441 
442  vComm->AllReduce(meshVertId, LibUtilities::ReduceSum);
443 
444  // When running in parallel, we need to ensure that the singular
445  // mesh vertex is communicated to any periodic vertices, otherwise
446  // the system may diverge.
447  if(m_systemSingular && checkIfSystemSingular)
448  {
449  // Firstly, we check that no other processors have this
450  // vertex. If they do, then we mark the vertex as also being
451  // Dirichlet.
452  if (maxIdx != p)
453  {
454  for (i = 0; i < locExpVector.size(); ++i)
455  {
456  for (j = 0; j < locExpVector[i]->GetNverts(); ++j)
457  {
458  if (locExpVector[i]->GetGeom()->GetVid(j) !=
459  meshVertId)
460  {
461  continue;
462  }
463 
464  if (graph[0].count(meshVertId) == 0)
465  {
466  graph[0][meshVertId] =
467  graphVertId++;
468  }
469  }
470  }
471  }
472 
473  // In the case that meshVertId is periodic with other vertices,
474  // this process and all other processes need to make sure that
475  // the periodic vertices are also marked as Dirichlet.
476  int gId;
477 
478  // At least one process (maxBCidx) will have already associated
479  // a graphVertId with meshVertId. Others won't even have any of
480  // the vertices. The logic below is designed to handle both
481  // cases.
482  if (graph[0].count(meshVertId) == 0)
483  {
484  gId = -1;
485  }
486  else
487  {
488  gId = graph[0][meshVertId];
489  }
490 
491  for (pIt = periodicVerts.begin();
492  pIt != periodicVerts.end(); ++pIt)
493  {
494  // Either the vertex is local to this processor (in which
495  // case it will be in the pIt->first position) or else
496  // meshVertId might be contained within another processor's
497  // vertex list. The if statement below covers both cases. If
498  // we find it, set as Dirichlet with the vertex id gId.
499  if (pIt->first == meshVertId)
500  {
501  gId = gId < 0 ? graphVertId++ : gId;
502  graph[0][meshVertId] = gId;
503 
504  for (i = 0; i < pIt->second.size(); ++i)
505  {
506  if (pIt->second[i].isLocal)
507  {
508  graph[0][pIt->second[i].id] = graph[0][meshVertId];
509  }
510  }
511  }
512  else
513  {
514  bool found = false;
515  for (i = 0; i < pIt->second.size(); ++i)
516  {
517  if (pIt->second[i].id == meshVertId)
518  {
519  found = true;
520  break;
521  }
522  }
523 
524  if (found)
525  {
526  gId = gId < 0 ? graphVertId++ : gId;
527  graph[0][pIt->first] = gId;
528 
529  for (i = 0; i < pIt->second.size(); ++i)
530  {
531  if (pIt->second[i].isLocal)
532  {
533  graph[0][pIt->second[i].id] = graph[0][pIt->first];
534  }
535  }
536  }
537  }
538  }
539  }
540 
541  // Add extra dirichlet boundary conditions to count.
542  m_numLocalDirBndCoeffs += nExtraDirichlet;
543  firstNonDirGraphVertId = graphVertId;
544 
545  typedef boost::adjacency_list<
546  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
547  BoostGraph boostGraphObj;
548 
549  vector<map<int,int> > tempGraph(3);
550  map<int, int> vwgts_map;
551  Array<OneD, int> localVerts;
552  Array<OneD, int> localEdges;
553  Array<OneD, int> localFaces;
554 
555  int tempGraphVertId = 0;
556  int localVertOffset = 0;
557  int localEdgeOffset = 0;
558  int localFaceOffset = 0;
559  int nTotalVerts = 0;
560  int nTotalEdges = 0;
561  int nTotalFaces = 0;
562  int nVerts;
563  int nEdges;
564  int nFaces;
565  int vertCnt;
566  int edgeCnt;
567  int faceCnt;
568 
569  m_numNonDirVertexModes = 0;
570  m_numNonDirEdges = 0;
571  m_numNonDirFaces = 0;
572  m_numNonDirFaceModes = 0;
573  m_numNonDirFaceModes = 0;
574  m_numLocalBndCoeffs = 0;
575 
576  map<int,int> EdgeSize;
577  map<int,int> FaceSize;
578 
579  /// - Count verts, edges, face and add up edges and face sizes
580  for(i = 0; i < locExpVector.size(); ++i)
581  {
582  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
583  nTotalVerts += exp->GetNverts();
584  nTotalEdges += exp->GetNedges();
585  nTotalFaces += exp->GetNfaces();
586 
587  nEdges = exp->GetNedges();
588  for(j = 0; j < nEdges; ++j)
589  {
590  meshEdgeId = exp->GetGeom()->GetEid(j);
591  if (EdgeSize.count(meshEdgeId) > 0)
592  {
593  EdgeSize[meshEdgeId] =
594  min(EdgeSize[meshEdgeId],
595  exp->GetEdgeNcoeffs(j) - 2);
596  }
597  else
598  {
599  EdgeSize[meshEdgeId] = exp->GetEdgeNcoeffs(j) - 2;
600  }
601  }
602 
603  nFaces = exp->GetNfaces();
604  faceCnt = 0;
605  for(j = 0; j < nFaces; ++j)
606  {
607  meshFaceId = exp->GetGeom()->GetFid(j);
608  if (FaceSize.count(meshFaceId) > 0)
609  {
610  FaceSize[meshFaceId] =
611  min(FaceSize[meshFaceId],
612  exp->GetFaceIntNcoeffs(j));
613  }
614  else
615  {
616  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
617  }
618  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
619  }
620  }
621 
622  /// - Periodic vertices
623  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
624  {
625  meshVertId = pIt->first;
626 
627  // This periodic vertex is joined to a Dirichlet condition.
628  if (graph[0].count(pIt->first) != 0)
629  {
630  for (i = 0; i < pIt->second.size(); ++i)
631  {
632  meshVertId2 = pIt->second[i].id;
633  if (graph[0].count(meshVertId2) == 0 &&
634  pIt->second[i].isLocal)
635  {
636  graph[0][meshVertId2] =
637  graph[0][meshVertId];
638  }
639  }
640  continue;
641  }
642 
643  // One of the attached vertices is Dirichlet.
644  bool isDirichlet = false;
645  for (i = 0; i < pIt->second.size(); ++i)
646  {
647  if (!pIt->second[i].isLocal)
648  {
649  continue;
650  }
651 
652  meshVertId2 = pIt->second[i].id;
653  if (graph[0].count(meshVertId2) > 0)
654  {
655  isDirichlet = true;
656  break;
657  }
658  }
659 
660  if (isDirichlet)
661  {
662  graph[0][meshVertId] =
663  graph[0][pIt->second[i].id];
664 
665  for (j = 0; j < pIt->second.size(); ++j)
666  {
667  meshVertId2 = pIt->second[i].id;
668  if (j == i || !pIt->second[j].isLocal ||
669  graph[0].count(meshVertId2) > 0)
670  {
671  continue;
672  }
673 
674  graph[0][meshVertId2] =
675  graph[0][pIt->second[i].id];
676  }
677 
678  continue;
679  }
680 
681  // Otherwise, see if a vertex ID has already been set.
682  for (i = 0; i < pIt->second.size(); ++i)
683  {
684  if (!pIt->second[i].isLocal)
685  {
686  continue;
687  }
688 
689  if (tempGraph[0].count(pIt->second[i].id) > 0)
690  {
691  break;
692  }
693  }
694 
695  if (i == pIt->second.size())
696  {
697  boost::add_vertex(boostGraphObj);
698  tempGraph[0][meshVertId] = tempGraphVertId++;
699  m_numNonDirVertexModes++;
700  }
701  else
702  {
703  tempGraph[0][meshVertId] = tempGraph[0][pIt->second[i].id];
704  }
705  }
706 
707  // Store the temporary graph vertex id's of all element edges and
708  // vertices in these 3 arrays below
709  localVerts = Array<OneD, int>(nTotalVerts,-1);
710  localEdges = Array<OneD, int>(nTotalEdges,-1);
711  localFaces = Array<OneD, int>(nTotalFaces,-1);
712 
713  // Set up vertex numbering
714  for(i = 0; i < locExpVector.size(); ++i)
715  {
716  exp = locExpVector[i];
717  vertCnt = 0;
718  nVerts = exp->GetNverts();
719  for(j = 0; j < nVerts; ++j)
720  {
721  meshVertId = exp->GetGeom()->GetVid(j);
722  if(graph[0].count(meshVertId) == 0)
723  {
724  if(tempGraph[0].count(meshVertId) == 0)
725  {
726  boost::add_vertex(boostGraphObj);
727  tempGraph[0][meshVertId] = tempGraphVertId++;
728  m_numNonDirVertexModes+=1;
729  }
730  localVerts[localVertOffset+vertCnt++] = tempGraph[0][meshVertId];
731  vwgts_map[ tempGraph[0][meshVertId] ] = 1;
732  }
733  }
734 
735  localVertOffset+=nVerts;
736  }
737 
738  /// - Periodic edges
739  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
740  {
741  meshEdgeId = pIt->first;
742 
743  // This periodic edge is joined to a Dirichlet condition.
744  if (graph[1].count(pIt->first) != 0)
745  {
746  for (i = 0; i < pIt->second.size(); ++i)
747  {
748  meshEdgeId2 = pIt->second[i].id;
749  if (graph[1].count(meshEdgeId2) == 0 &&
750  pIt->second[i].isLocal)
751  {
752  graph[1][meshEdgeId2] =
753  graph[1][meshEdgeId];
754  }
755  }
756  continue;
757  }
758 
759  // One of the attached edges is Dirichlet.
760  bool isDirichlet = false;
761  for (i = 0; i < pIt->second.size(); ++i)
762  {
763  if (!pIt->second[i].isLocal)
764  {
765  continue;
766  }
767 
768  meshEdgeId2 = pIt->second[i].id;
769  if (graph[1].count(meshEdgeId2) > 0)
770  {
771  isDirichlet = true;
772  break;
773  }
774  }
775 
776  if (isDirichlet)
777  {
778  graph[1][meshEdgeId] =
779  graph[1][pIt->second[i].id];
780 
781  for (j = 0; j < pIt->second.size(); ++j)
782  {
783  meshEdgeId2 = pIt->second[i].id;
784  if (j == i || !pIt->second[j].isLocal ||
785  graph[1].count(meshEdgeId2) > 0)
786  {
787  continue;
788  }
789 
790  graph[1][meshEdgeId2] =
791  graph[1][pIt->second[i].id];
792  }
793 
794  continue;
795  }
796 
797  // Otherwise, see if a edge ID has already been set.
798  for (i = 0; i < pIt->second.size(); ++i)
799  {
800  if (!pIt->second[i].isLocal)
801  {
802  continue;
803  }
804 
805  if (tempGraph[1].count(pIt->second[i].id) > 0)
806  {
807  break;
808  }
809  }
810 
811  if (i == pIt->second.size())
812  {
813  boost::add_vertex(boostGraphObj);
814  tempGraph[1][meshEdgeId] = tempGraphVertId++;
815  m_numNonDirEdgeModes += EdgeSize[meshEdgeId];
816  m_numNonDirEdges++;
817  }
818  else
819  {
820  tempGraph[1][meshEdgeId] = tempGraph[1][pIt->second[i].id];
821  }
822  }
823 
824  int nEdgeIntCoeffs, nFaceIntCoeffs;
825 
826  // Set up edge numbering
827  for(i = 0; i < locExpVector.size(); ++i)
828  {
829  exp = locExpVector[i];
830  edgeCnt = 0;
831  nEdges = exp->GetNedges();
832 
833  for(j = 0; j < nEdges; ++j)
834  {
835  nEdgeIntCoeffs = exp->GetEdgeNcoeffs(j) - 2;
836  meshEdgeId = exp->GetGeom()->GetEid(j);
837  if(graph[1].count(meshEdgeId) == 0)
838  {
839  if(tempGraph[1].count(meshEdgeId) == 0)
840  {
841  boost::add_vertex(boostGraphObj);
842  tempGraph[1][meshEdgeId] = tempGraphVertId++;
843  m_numNonDirEdgeModes+=nEdgeIntCoeffs;
844 
845  m_numNonDirEdges++;
846  }
847  localEdges[localEdgeOffset+edgeCnt++] = tempGraph[1][meshEdgeId];
848  vwgts_map[ tempGraph[1][meshEdgeId] ] = nEdgeIntCoeffs;
849  }
850  }
851 
852  localEdgeOffset+=nEdges;
853  }
854 
855  /// - Periodic faces
856  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
857  {
858  if (!pIt->second[0].isLocal)
859  {
860  // The face mapped to is on another process.
861  meshFaceId = pIt->first;
862  ASSERTL0(graph[2].count(meshFaceId) == 0,
863  "This periodic boundary edge has been specified before");
864  boost::add_vertex(boostGraphObj);
865  tempGraph[2][meshFaceId] = tempGraphVertId++;
866  nFaceIntCoeffs = FaceSize[meshFaceId];
867  m_numNonDirFaceModes+=nFaceIntCoeffs;
868  m_numNonDirFaces++;
869  }
870  else if (pIt->first < pIt->second[0].id)
871  {
872  ASSERTL0(graph[2].count(pIt->first) == 0,
873  "This periodic boundary face has been specified before");
874  ASSERTL0(graph[2].count(pIt->second[0].id) == 0,
875  "This periodic boundary face has been specified before");
876 
877  boost::add_vertex(boostGraphObj);
878  tempGraph[2][pIt->first] = tempGraphVertId;
879  tempGraph[2][pIt->second[0].id] = tempGraphVertId++;
880  nFaceIntCoeffs = FaceSize[pIt->first];
881  m_numNonDirFaceModes+=nFaceIntCoeffs;
882  m_numNonDirFaces++;
883  }
884  }
885 
886  // setup face numbering
887  for(i = 0; i < locExpVector.size(); ++i)
888  {
889  exp = locExpVector[i];
890  nFaces = exp->GetNfaces();
891  faceCnt = 0;
892  for(j = 0; j < nFaces; ++j)
893  {
894  nFaceIntCoeffs = exp->GetFaceIntNcoeffs(j);
895  meshFaceId = exp->GetGeom()->GetFid(j);
896  if(graph[2].count(meshFaceId) == 0)
897  {
898  if(tempGraph[2].count(meshFaceId) == 0)
899  {
900  boost::add_vertex(boostGraphObj);
901  tempGraph[2][meshFaceId] = tempGraphVertId++;
902  m_numNonDirFaceModes+=nFaceIntCoeffs;
903 
904  m_numNonDirFaces++;
905  }
906  localFaces[localFaceOffset+faceCnt++] = tempGraph[2][meshFaceId];
907  vwgts_map[ tempGraph[2][meshFaceId] ] = nFaceIntCoeffs;
908  }
909  }
910  m_numLocalBndCoeffs += exp->NumBndryCoeffs();
911 
912  localFaceOffset+=nFaces;
913  }
914 
915  localVertOffset=0;
916  localEdgeOffset=0;
917  localFaceOffset=0;
918  for(i = 0; i < locExpVector.size(); ++i)
919  {
920  exp = locExpVector[i];
921  nVerts = exp->GetNverts();
922  nEdges = exp->GetNedges();
923  nFaces = exp->GetNfaces();
924 
925  // Now loop over all local faces, edges and vertices of this
926  // element and define that all other faces, edges and verices of
927  // this element are adjacent to them.
928 
929  // Vertices
930  for(j = 0; j < nVerts; j++)
931  {
932  if(localVerts[j+localVertOffset]==-1)
933  {
934  break;
935  }
936  // associate to other vertices
937  for(k = 0; k < nVerts; k++)
938  {
939  if(localVerts[k+localVertOffset]==-1)
940  {
941  break;
942  }
943  if(k!=j)
944  {
945  boost::add_edge( (size_t) localVerts[j+localVertOffset],
946  (size_t) localVerts[k+localVertOffset],boostGraphObj);
947  }
948  }
949  // associate to other edges
950  for(k = 0; k < nEdges; k++)
951  {
952  if(localEdges[k+localEdgeOffset]==-1)
953  {
954  break;
955  }
956  boost::add_edge( (size_t) localVerts[j+localVertOffset],
957  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
958  }
959  // associate to other faces
960  for(k = 0; k < nFaces; k++)
961  {
962  if(localFaces[k+localFaceOffset]==-1)
963  {
964  break;
965  }
966  boost::add_edge( (size_t) localVerts[j+localVertOffset],
967  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
968  }
969  }
970 
971  // Edges
972  for(j = 0; j < nEdges; j++)
973  {
974  if(localEdges[j+localEdgeOffset]==-1)
975  {
976  break;
977  }
978  // Associate to other edges
979  for(k = 0; k < nEdges; k++)
980  {
981  if(localEdges[k+localEdgeOffset]==-1)
982  {
983  break;
984  }
985  if(k!=j)
986  {
987  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
988  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
989  }
990  }
991  // Associate to vertices
992  for(k = 0; k < nVerts; k++)
993  {
994  if(localVerts[k+localVertOffset]==-1)
995  {
996  break;
997  }
998  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
999  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1000  }
1001  // Associate to faces
1002  for(k = 0; k < nFaces; k++)
1003  {
1004  if(localFaces[k+localFaceOffset]==-1)
1005  {
1006  break;
1007  }
1008  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1009  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1010  }
1011  }
1012 
1013  // Faces
1014  for(j = 0; j < nFaces; j++)
1015  {
1016  if(localFaces[j+localFaceOffset]==-1)
1017  {
1018  break;
1019  }
1020  // Associate to other faces
1021  for(k = 0; k < nFaces; k++)
1022  {
1023  if(localFaces[k+localFaceOffset]==-1)
1024  {
1025  break;
1026  }
1027  if(k!=j)
1028  {
1029  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1030  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1031  }
1032  }
1033  // Associate to vertices
1034  for(k = 0; k < nVerts; k++)
1035  {
1036  if(localVerts[k+localVertOffset]==-1)
1037  {
1038  break;
1039  }
1040  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1041  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1042  }
1043  // Associate to edges
1044  for(k = 0; k < nEdges; k++)
1045  {
1046  if(localEdges[k+localEdgeOffset]==-1)
1047  {
1048  break;
1049  }
1050  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1051  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1052  }
1053  }
1054 
1055  localVertOffset+=nVerts;
1056  localEdgeOffset+=nEdges;
1057  localFaceOffset+=nFaces;
1058  }
1059 
1060  // Container to store vertices of the graph which correspond to
1061  // degrees of freedom along the boundary and periodic BCs.
1062  set<int> partVerts;
1063 
1064  if (m_solnType == eIterativeMultiLevelStaticCond ||
1065  m_solnType == eXxtMultiLevelStaticCond)
1066  {
1067  vector<long> procVerts, procEdges, procFaces;
1068  set <int> foundVerts, foundEdges, foundFaces;
1069 
1070  // Loop over element and construct the procVerts and procEdges
1071  // vectors, which store the geometry IDs of mesh vertices and
1072  // edges respectively which are local to this process.
1073  for(i = cnt = 0; i < locExpVector.size(); ++i)
1074  {
1075  int elmtid = locExp.GetOffset_Elmt_Id(i);
1076  exp = locExpVector[elmtid];
1077  for (j = 0; j < exp->GetNverts(); ++j)
1078  {
1079  int vid = exp->GetGeom()->GetVid(j)+1;
1080  if (foundVerts.count(vid) == 0)
1081  {
1082  procVerts.push_back(vid);
1083  foundVerts.insert(vid);
1084  }
1085  }
1086 
1087  for (j = 0; j < exp->GetNedges(); ++j)
1088  {
1089  int eid = exp->GetGeom()->GetEid(j)+1;
1090 
1091  if (foundEdges.count(eid) == 0)
1092  {
1093  procEdges.push_back(eid);
1094  foundEdges.insert(eid);
1095  }
1096  }
1097 
1098  for (j = 0; j < exp->GetNfaces(); ++j)
1099  {
1100  int fid = exp->GetGeom()->GetFid(j)+1;
1101 
1102  if (foundFaces.count(fid) == 0)
1103  {
1104  procFaces.push_back(fid);
1105  foundFaces.insert(fid);
1106  }
1107  }
1108  }
1109 
1110  int unique_verts = foundVerts.size();
1111  int unique_edges = foundEdges.size();
1112  int unique_faces = foundFaces.size();
1113 
1114  // Now construct temporary GS objects. These will be used to
1115  // populate the arrays tmp3 and tmp4 with the multiplicity of
1116  // the vertices and edges respectively to identify those
1117  // vertices and edges which are located on partition boundary.
1118  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1119  Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm);
1120  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1121  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1122  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1123  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1124  Gs::Finalise(tmp1);
1125 
1126  if (unique_edges > 0)
1127  {
1128  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1129  Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
1130  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1131  Gs::Finalise(tmp2);
1132  }
1133 
1134  if (unique_faces > 0)
1135  {
1136  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1137  Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
1138  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1139  Gs::Finalise(tmp3);
1140  }
1141 
1142  // Finally, fill the partVerts set with all non-Dirichlet
1143  // vertices which lie on a partition boundary.
1144  for (i = 0; i < unique_verts; ++i)
1145  {
1146  if (tmp4[i] > 1.0)
1147  {
1148  if (graph[0].count(procVerts[i]-1) == 0)
1149  {
1150  partVerts.insert(tempGraph[0][procVerts[i]-1]);
1151  }
1152  }
1153  }
1154 
1155  for (i = 0; i < unique_edges; ++i)
1156  {
1157  if (tmp5[i] > 1.0)
1158  {
1159  if (graph[1].count(procEdges[i]-1) == 0)
1160  {
1161  partVerts.insert(tempGraph[1][procEdges[i]-1]);
1162  }
1163  }
1164  }
1165 
1166  for (i = 0; i < unique_faces; ++i)
1167  {
1168  if (tmp6[i] > 1.0)
1169  {
1170  if (graph[2].count(procFaces[i]-1) == 0)
1171  {
1172  partVerts.insert(tempGraph[2][procFaces[i]-1]);
1173  }
1174  }
1175  }
1176 
1177  // Now fill with all vertices on periodic BCs
1178  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
1179  {
1180  if (graph[0].count(pIt->first) == 0)
1181  {
1182  partVerts.insert(tempGraph[0][pIt->first]);
1183  }
1184  }
1185  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1186  {
1187  if (graph[1].count(pIt->first) == 0)
1188  {
1189  partVerts.insert(tempGraph[1][pIt->first]);
1190  }
1191  }
1192  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1193  {
1194  if (graph[2].count(pIt->first) == 0)
1195  {
1196  partVerts.insert(tempGraph[2][pIt->first]);
1197  }
1198  }
1199  }
1200 
1201  int nGraphVerts = tempGraphVertId;
1202  Array<OneD, int> perm(nGraphVerts);
1203  Array<OneD, int> iperm(nGraphVerts);
1204  Array<OneD, int> vwgts(nGraphVerts);
1205  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1206  for(i = 0; i < nGraphVerts; ++i)
1207  {
1208  vwgts[i] = vwgts_map[i];
1209  }
1210 
1211  if(nGraphVerts)
1212  {
1213  switch(m_solnType)
1214  {
1215  case eDirectFullMatrix:
1216  case eIterativeFull:
1217  case eIterativeStaticCond:
1218  case ePETScStaticCond:
1219  case ePETScFullMatrix:
1220  case eXxtFullMatrix:
1221  case eXxtStaticCond:
1222  {
1223  NoReordering(boostGraphObj,perm,iperm);
1224  break;
1225  }
1226 
1227  case eDirectStaticCond:
1228  {
1229  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1230  break;
1231  }
1232 
1233  case ePETScMultiLevelStaticCond:
1234  case eDirectMultiLevelStaticCond:
1235  case eIterativeMultiLevelStaticCond:
1236  case eXxtMultiLevelStaticCond:
1237  {
1238  MultiLevelBisectionReordering(
1239  boostGraphObj, perm, iperm, bottomUpGraph,
1240  partVerts, mdswitch);
1241  break;
1242  }
1243  default:
1244  {
1245  ASSERTL0(false,
1246  "Unrecognised solution type: " + std::string(
1247  GlobalSysSolnTypeMap[m_solnType]));
1248  }
1249  }
1250  }
1251 
1252  // For parallel multi-level static condensation determine the lowest
1253  // static condensation level amongst processors.
1254  if (m_solnType == eIterativeMultiLevelStaticCond ||
1255  m_solnType == eXxtMultiLevelStaticCond)
1256  {
1257  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1258  vComm->AllReduce(m_lowestStaticCondLevel,
1259  LibUtilities::ReduceMax);
1260  }
1261  else
1262  {
1263  m_lowestStaticCondLevel = 0;
1264  }
1265 
1266  /**
1267  * STEP 4: Fill the #graph[0] and
1268  * #graph[1] with the optimal ordering from boost.
1269  */
1270  map<int, int>::iterator mapIt;
1271  for(mapIt = tempGraph[0].begin(); mapIt != tempGraph[0].end(); mapIt++)
1272  {
1273  graph[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1274  }
1275  for(mapIt = tempGraph[1].begin(); mapIt != tempGraph[1].end(); mapIt++)
1276  {
1277  graph[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1278  }
1279  for(mapIt = tempGraph[2].begin(); mapIt != tempGraph[2].end(); mapIt++)
1280  {
1281  graph[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1282  }
1283 
1284  return nGraphVerts;
1285  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MultiRegions::AssemblyMap::m_systemSingular
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:320
Nektar::MultiRegions::MultiLevelBisectionReordering
void MultiLevelBisectionReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm, BottomUpSubStructuredGraphSharedPtr &substructgraph, std::set< int > partVerts, int mdswitch)
Definition: SubStructuredGraph.cpp:862
Nektar::MultiRegions::AssemblyMap::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
Gs::Gather
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:224
Gs::Finalise
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:206
Nektar::MultiRegions::AssemblyMapCG::m_extraDirEdges
Array< OneD, int > m_extraDirEdges
Extra dirichlet edges in parallel.
Definition: AssemblyMapCG.h:141
Gs::Init
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirVertexModes
int m_numNonDirVertexModes
Number of non Dirichlet vertex modes.
Definition: AssemblyMapCG.h:125
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdges
int m_numNonDirEdges
Number of Dirichlet edges.
Definition: AssemblyMapCG.h:135
Nektar::MultiRegions::CuthillMckeeReordering
void CuthillMckeeReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
Definition: SubStructuredGraph.cpp:839
Nektar::MultiRegions::AssemblyMapCG::m_numLocalBndCondCoeffs
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
Definition: AssemblyMapCG.h:139
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:741
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Nektar::MultiRegions::GlobalSysSolnTypeMap
const char *const GlobalSysSolnTypeMap[]
Definition: MultiRegions.hpp:96
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaceModes
int m_numNonDirFaceModes
Number of non Dirichlet face modes.
Definition: AssemblyMapCG.h:129
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::MultiRegions::AssemblyMap::m_solnType
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:363
Nektar::MultiRegions::AssemblyMapCG::m_numDirFaces
int m_numDirFaces
Number of Dirichlet faces.
Definition: AssemblyMapCG.h:133
Nektar::LocalRegions::ExpansionSharedPtr
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::MultiRegions::AssemblyMap::m_lowestStaticCondLevel
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:397
Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:316
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::AssemblyMapCG::m_numDirEdges
int m_numDirEdges
Number of Dirichlet edges.
Definition: AssemblyMapCG.h:131
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdgeModes
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.
Definition: AssemblyMapCG.h:127
Gs::gs_data
Definition: GsLib.hpp:117
Nektar::MultiRegions::AssemblyMap::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:303
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaces
int m_numNonDirFaces
Number of Dirichlet faces.
Definition: AssemblyMapCG.h:137
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
Nektar::MultiRegions::NoReordering
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
Definition: SubStructuredGraph.cpp:1074
std::map<int, std::vector<ExtraDirDof> >& Nektar::MultiRegions::AssemblyMapCG::GetExtraDirDofs ( )
inline

Definition at line 108 of file AssemblyMapCG.h.

References m_extraDirDofs.

109  {
110  return m_extraDirDofs;
111  }
Nektar::MultiRegions::AssemblyMapCG::m_extraDirDofs
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
Definition: AssemblyMapCG.h:148
void Nektar::MultiRegions::AssemblyMapCG::SetUpUniversalC0ContMap ( const ExpList &  locExp,
const PeriodicMap &  perVerts = NullPeriodicMap,
const PeriodicMap &  perEdges = NullPeriodicMap,
const PeriodicMap &  perFaces = NullPeriodicMap 
)
protected

Sets up the global to universal mapping of degrees of freedom across processors.

Definition at line 1995 of file AssemblyMapCG.cpp.

References Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Gs::Init(), Nektar::MultiRegions::AssemblyMap::m_bndGsh, Nektar::MultiRegions::AssemblyMap::m_comm, Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMap, Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMapUnique, m_globalToUniversalMap, m_globalToUniversalMapUnique, Nektar::MultiRegions::AssemblyMap::m_gsh, m_localToGlobalMap, Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Gs::Unique(), and Vmath::Zero().

Referenced by AssemblyMapCG().

2000  {
2001  LocalRegions::ExpansionSharedPtr exp;
2002  int nVert = 0;
2003  int nEdge = 0;
2004  int nFace = 0;
2005  int maxEdgeDof = 0;
2006  int maxFaceDof = 0;
2007  int maxIntDof = 0;
2008  int dof = 0;
2009  int cnt;
2010  int i,j,k;
2011  int meshVertId;
2012  int meshEdgeId;
2013  int meshFaceId;
2014  int elementId;
2015  int vGlobalId;
2016  int maxBndGlobalId = 0;
2017  StdRegions::Orientation edgeOrient;
2018  StdRegions::Orientation faceOrient;
2019  Array<OneD, unsigned int> edgeInteriorMap;
2020  Array<OneD, int> edgeInteriorSign;
2021  Array<OneD, unsigned int> faceInteriorMap;
2022  Array<OneD, int> faceInteriorSign;
2023  Array<OneD, unsigned int> interiorMap;
2024  PeriodicMap::const_iterator pIt;
2025 
2026  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
2027  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
2028 
2029  m_globalToUniversalMap = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
2030  m_globalToUniversalMapUnique = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
2031  m_globalToUniversalBndMap = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
2032  m_globalToUniversalBndMapUnique = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
2033 
2034  // Loop over all the elements in the domain to gather mesh data
2035  for(i = 0; i < locExpVector.size(); ++i)
2036  {
2037  exp = locExpVector[i];
2038  nVert += exp->GetNverts();
2039  nEdge += exp->GetNedges();
2040  nFace += exp->GetNfaces();
2041  // Loop over all edges (and vertices) of element i
2042  for(j = 0; j < exp->GetNedges(); ++j)
2043  {
2044  dof = exp->GetEdgeNcoeffs(j)-2;
2045  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2046  }
2047  for(j = 0; j < exp->GetNfaces(); ++j)
2048  {
2049  dof = exp->GetFaceIntNcoeffs(j);
2050  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2051  }
2052  exp->GetInteriorMap(interiorMap);
2053  dof = interiorMap.num_elements();
2054  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2055  }
2056 
2057  // Tell other processes about how many dof we have
2058  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2059  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2060  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2061  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2062  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2063  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2064 
2065  // Assemble global to universal mapping for this process
2066  for(i = 0; i < locExpVector.size(); ++i)
2067  {
2068  exp = locExpVector[i];
2069  cnt = locExp.GetCoeff_Offset(i);
2070 
2071  // Loop over all vertices of element i
2072  for(j = 0; j < exp->GetNverts(); ++j)
2073  {
2074  meshVertId = exp->GetGeom()->GetVid(j);
2075  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2076 
2077  pIt = perVerts.find(meshVertId);
2078  if (pIt != perVerts.end())
2079  {
2080  for (k = 0; k < pIt->second.size(); ++k)
2081  {
2082  meshVertId = min(meshVertId, pIt->second[k].id);
2083  }
2084  }
2085 
2086  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2087  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2088  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2089  }
2090 
2091  // Loop over all edges of element i
2092  for(j = 0; j < exp->GetNedges(); ++j)
2093  {
2094  meshEdgeId = exp->GetGeom()->GetEid(j);
2095  pIt = perEdges.find(meshEdgeId);
2096  edgeOrient = exp->GetGeom()->GetEorient(j);
2097 
2098  if (pIt != perEdges.end())
2099  {
2100  pair<int, StdRegions::Orientation> idOrient =
2101  DeterminePeriodicEdgeOrientId(
2102  meshEdgeId, edgeOrient, pIt->second);
2103  meshEdgeId = idOrient.first;
2104  edgeOrient = idOrient.second;
2105  }
2106 
2107  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
2108  dof = exp->GetEdgeNcoeffs(j)-2;
2109 
2110  // Set the global DOF's for the interior modes of edge j
2111  // run backwards because of variable P case "ghost" modes
2112  for(k = dof-1; k >= 0; --k)
2113  {
2114  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2115  m_globalToUniversalMap[vGlobalId]
2116  = nVert + meshEdgeId * maxEdgeDof + k + 1;
2117  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2118  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2119  }
2120  }
2121 
2122  // Loop over all faces of element i
2123  for(j = 0; j < exp->GetNfaces(); ++j)
2124  {
2125  faceOrient = exp->GetGeom()->GetForient(j);
2126 
2127  meshFaceId = exp->GetGeom()->GetFid(j);
2128 
2129  pIt = perFaces.find(meshFaceId);
2130  if (pIt != perFaces.end())
2131  {
2132  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2133  {
2134  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
2135  }
2136  meshFaceId = min(meshFaceId, pIt->second[0].id);
2137  }
2138 
2139 
2140  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2141  dof = exp->GetFaceIntNcoeffs(j);
2142 
2143 
2144  for(k = dof-1; k >= 0; --k)
2145  {
2146  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2147  m_globalToUniversalMap[vGlobalId]
2148  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2149  + k + 1;
2150  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2151 
2152  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2153  }
2154  }
2155 
2156  // Add interior DOFs to complete universal numbering
2157  exp->GetInteriorMap(interiorMap);
2158  dof = interiorMap.num_elements();
2159  elementId = (exp->GetGeom())->GetGlobalID();
2160  for (k = 0; k < dof; ++k)
2161  {
2162  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2163  m_globalToUniversalMap[vGlobalId]
2164  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2165  }
2166  }
2167 
2168  // Set up the GSLib universal assemble mapping
2169  // Internal DOF do not participate in any data
2170  // exchange, so we keep these set to the special GSLib id=0 so
2171  // they are ignored.
2172  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2173  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2174  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2175  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2176  {
2177  tmp[i] = m_globalToUniversalMap[i];
2178  }
2179 
2180  m_gsh = Gs::Init(tmp, vCommRow);
2181  m_bndGsh = Gs::Init(tmp2, vCommRow);
2182  Gs::Unique(tmp, vCommRow);
2183  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2184  {
2185  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2186  }
2187  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2188  {
2189  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2190  }
2191  }
Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
Nektar::MultiRegions::AssemblyMap::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMapUnique
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMapCG.h:123
Gs::Init
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:121
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::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::MultiRegions::DeterminePeriodicFaceOrient
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:1952
Gs::Unique
static void Unique(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Updates pId to negate all-but-one references to each universal ID.
Definition: GsLib.hpp:184
Nektar::LocalRegions::ExpansionSharedPtr
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMap
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:358
Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMapUnique
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:360
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:1904
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Nektar::MultiRegions::AssemblyMapCG::v_Assemble ( const Array< OneD, const NekDouble > &  loc,
Array< OneD, NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2495 of file AssemblyMapCG.cpp.

References Vmath::Assmb(), m_localToGlobalMap, m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::UniversalAssemble(), and Vmath::Zero().

2498  {
2499  Array<OneD, const NekDouble> local;
2500  if(global.data() == loc.data())
2501  {
2502  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2503  }
2504  else
2505  {
2506  local = loc; // create reference
2507  }
2508  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2509 
2510  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2511 
2512  if(m_signChange)
2513  {
2514  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2515  }
2516  else
2517  {
2518  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2519  }
2520  UniversalAssemble(global);
2521  }
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
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::UniversalAssemble
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:764
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:695
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:117
Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Nektar::MultiRegions::AssemblyMapCG::v_Assemble ( const NekVector< NekDouble > &  loc,
NekVector< NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2523 of file AssemblyMapCG.cpp.

References Nektar::MultiRegions::AssemblyMap::Assemble(), and Nektar::NekVector< DataType >::GetPtr().

2526  {
2527  Assemble(loc.GetPtr(),global.GetPtr());
2528  }
Nektar::MultiRegions::AssemblyMap::Assemble
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Definition: AssemblyMap.cpp:750
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetExtraDirEdges ( )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2592 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2593  {
2594  return m_extraDirEdges;
2595  }
Nektar::MultiRegions::AssemblyMapCG::m_extraDirEdges
Array< OneD, int > m_extraDirEdges
Extra dirichlet edges in parallel.
Definition: AssemblyMapCG.h:141
int Nektar::MultiRegions::AssemblyMapCG::v_GetFullSystemBandWidth ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2552 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2553  {
2554  return m_fullSystemBandWidth;
2555  }
Nektar::MultiRegions::AssemblyMapCG::m_fullSystemBandWidth
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
Definition: AssemblyMapCG.h:119
int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2382 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2383  {
2384  return m_globalToUniversalMap[i];
2385  }
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:121
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2399 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2400  {
2401  return m_globalToUniversalMap;
2402  }
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:121
int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2387 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2388  {
2389  return m_globalToUniversalMapUnique[i];
2390  }
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMapUnique
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMapCG.h:123
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2405 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2406  {
2407  return m_globalToUniversalMapUnique;
2408  }
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMapUnique
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMapCG.h:123
int Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2377 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2378  {
2379  return m_localToGlobalMap[i];
2380  }
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:115
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2393 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2394  {
2395  return m_localToGlobalMap;
2396  }
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:115
NekDouble Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2410 of file AssemblyMapCG.cpp.

References m_localToGlobalSign, and Nektar::MultiRegions::AssemblyMap::m_signChange.

2412  {
2413  if(m_signChange)
2414  {
2415  return m_localToGlobalSign[i];
2416  }
2417  else
2418  {
2419  return 1.0;
2420  }
2421  }
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:117
const Array< OneD, NekDouble > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2423 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2424  {
2425  return m_localToGlobalSign;
2426  }
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:117
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2572 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2573  {
2574  return m_numDirEdges;
2575  }
Nektar::MultiRegions::AssemblyMapCG::m_numDirEdges
int m_numDirEdges
Number of Dirichlet edges.
Definition: AssemblyMapCG.h:131
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2577 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2578  {
2579  return m_numDirFaces;
2580  }
Nektar::MultiRegions::AssemblyMapCG::m_numDirFaces
int m_numDirFaces
Number of Dirichlet faces.
Definition: AssemblyMapCG.h:133
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2562 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2563  {
2564  return m_numNonDirEdgeModes;
2565  }
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdgeModes
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.
Definition: AssemblyMapCG.h:127
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2582 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2583  {
2584  return m_numNonDirEdges;
2585  }
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdges
int m_numNonDirEdges
Number of Dirichlet edges.
Definition: AssemblyMapCG.h:135
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2567 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2568  {
2569  return m_numNonDirFaceModes;
2570  }
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaceModes
int m_numNonDirFaceModes
Number of non Dirichlet face modes.
Definition: AssemblyMapCG.h:129
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2587 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2588  {
2589  return m_numNonDirFaces;
2590  }
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaces
int m_numNonDirFaces
Number of Dirichlet faces.
Definition: AssemblyMapCG.h:137
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2557 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2558  {
2559  return m_numNonDirVertexModes;
2560  }
Nektar::MultiRegions::AssemblyMapCG::m_numNonDirVertexModes
int m_numNonDirVertexModes
Number of non Dirichlet vertex modes.
Definition: AssemblyMapCG.h:125
void Nektar::MultiRegions::AssemblyMapCG::v_GlobalToLocal ( const Array< OneD, const NekDouble > &  global,
Array< OneD, NekDouble > &  loc 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2463 of file AssemblyMapCG.cpp.

References Vmath::Gathr(), m_localToGlobalMap, m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, and Nektar::MultiRegions::AssemblyMap::m_signChange.

2466  {
2467  Array<OneD, const NekDouble> glo;
2468  if(global.data() == loc.data())
2469  {
2470  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2471  }
2472  else
2473  {
2474  glo = global; // create reference
2475  }
2476 
2477 
2478  if(m_signChange)
2479  {
2480  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2481  }
2482  else
2483  {
2484  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2485  }
2486  }
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
Vmath::Gathr
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:115
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:117
void Nektar::MultiRegions::AssemblyMapCG::v_GlobalToLocal ( const NekVector< NekDouble > &  global,
NekVector< NekDouble > &  loc 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2488 of file AssemblyMapCG.cpp.

References Nektar::NekVector< DataType >::GetPtr(), and Nektar::MultiRegions::AssemblyMap::GlobalToLocal().

2491  {
2492  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2493  }
Nektar::MultiRegions::AssemblyMap::GlobalToLocal
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
Definition: AssemblyMap.cpp:736
AssemblyMapSharedPtr Nektar::MultiRegions::AssemblyMapCG::v_LinearSpaceMap ( const ExpList &  locexp,
GlobalSysSolnType  solnType 
)
protectedvirtual

Construct an AssemblyMapCG object which corresponds to the linear space of the current object.

This function is used to create a linear-space assembly map, which is then used in the linear space preconditioner in the conjugate gradient solve.

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2201 of file AssemblyMapCG.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::eNull, Nektar::MultiRegions::ExpList::GetExp(), Gs::Init(), Nektar::MultiRegions::AssemblyMap::m_comm, m_globalToUniversalMap, m_localToGlobalMap, Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_session, Gs::Unique(), and Vmath::Zero().

2203  {
2204  AssemblyMapCGSharedPtr returnval;
2205 
2206  int i, j;
2207  int nverts = 0;
2208  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2209  = locexp.GetExp();
2210  int nelmts = exp->size();
2211 
2212  // Get Default Map and turn off any searched values.
2213  returnval = MemoryManager<AssemblyMapCG>
2214  ::AllocateSharedPtr(m_session);
2215  returnval->m_solnType = solnType;
2216  returnval->m_preconType = eNull;
2217  returnval->m_maxStaticCondLevel = 0;
2218  returnval->m_signChange = false;
2219  returnval->m_comm = m_comm;
2220 
2221  // Count the number of vertices
2222  for (i = 0; i < nelmts; ++i)
2223  {
2224  nverts += (*exp)[i]->GetNverts();
2225  }
2226 
2227  returnval->m_numLocalCoeffs = nverts;
2228  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2229 
2230  // Store original global ids in this map
2231  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2232 
2233  int cnt = 0;
2234  int cnt1 = 0;
2235  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2236 
2237  // Set up local to global map;
2238  for (i = 0; i < nelmts; ++i)
2239  {
2240  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2241  {
2242  returnval->m_localToGlobalMap[cnt] =
2243  returnval->m_localToGlobalBndMap[cnt] =
2244  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2245  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2246 
2247  // Set up numLocalDirBndCoeffs
2248  if ((returnval->m_localToGlobalMap[cnt]) <
2249  m_numGlobalDirBndCoeffs)
2250  {
2251  returnval->m_numLocalDirBndCoeffs++;
2252  }
2253  cnt++;
2254  }
2255  cnt1 += (*exp)[i]->GetNcoeffs();
2256  }
2257 
2258  cnt = 0;
2259  // Reset global numbering and count number of dofs
2260  for (i = 0; i < m_numGlobalCoeffs; ++i)
2261  {
2262  if (GlobCoeffs[i] != -1)
2263  {
2264  GlobCoeffs[i] = cnt++;
2265  }
2266  }
2267 
2268  // Set up number of globalCoeffs;
2269  returnval->m_numGlobalCoeffs = cnt;
2270 
2271  // Set up number of global Dirichlet boundary coefficients
2272  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2273  {
2274  if (GlobCoeffs[i] != -1)
2275  {
2276  returnval->m_numGlobalDirBndCoeffs++;
2277  }
2278  }
2279 
2280  // Set up global to universal map
2281  if (m_globalToUniversalMap.num_elements())
2282  {
2283  LibUtilities::CommSharedPtr vCommRow
2284  = m_session->GetComm()->GetRowComm();
2285  int nglocoeffs = returnval->m_numGlobalCoeffs;
2286  returnval->m_globalToUniversalMap
2287  = Array<OneD, int> (nglocoeffs);
2288  returnval->m_globalToUniversalMapUnique
2289  = Array<OneD, int> (nglocoeffs);
2290 
2291  // Reset local to global map and setup universal map
2292  for (i = 0; i < nverts; ++i)
2293  {
2294  cnt = returnval->m_localToGlobalMap[i];
2295  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2296 
2297  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2298  m_globalToUniversalMap[cnt];
2299  }
2300 
2301  Nektar::Array<OneD, long> tmp(nglocoeffs);
2302  Vmath::Zero(nglocoeffs, tmp, 1);
2303  for (unsigned int i = 0; i < nglocoeffs; ++i)
2304  {
2305  tmp[i] = returnval->m_globalToUniversalMap[i];
2306  }
2307  returnval->m_gsh = Gs::Init(tmp, vCommRow);
2308  Gs::Unique(tmp, vCommRow);
2309  for (unsigned int i = 0; i < nglocoeffs; ++i)
2310  {
2311  returnval->m_globalToUniversalMapUnique[i]
2312  = (tmp[i] >= 0 ? 1 : 0);
2313  }
2314  }
2315  else // not sure this option is ever needed.
2316  {
2317  for (i = 0; i < nverts; ++i)
2318  {
2319  cnt = returnval->m_localToGlobalMap[i];
2320  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2321  }
2322  }
2323 
2324  return returnval;
2325  }
Nektar::MultiRegions::AssemblyMap::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Gs::Init
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:121
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:115
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
Gs::Unique
static void Unique(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Updates pId to negate all-but-one references to each universal ID.
Definition: GsLib.hpp:184
Nektar::MultiRegions::eNull
No Solution type specified.
Definition: MultiRegions.hpp:123
Nektar::MultiRegions::AssemblyMap::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:303
Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::MultiRegions::AssemblyMapCGSharedPtr
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:52
void Nektar::MultiRegions::AssemblyMapCG::v_LocalToGlobal ( const Array< OneD, const NekDouble > &  loc,
Array< OneD, NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2428 of file AssemblyMapCG.cpp.

References Gs::Gather(), Gs::gs_max, Nektar::MultiRegions::AssemblyMap::m_gsh, m_localToGlobalMap, m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_signChange, and Vmath::Scatr().

2431  {
2432  Array<OneD, const NekDouble> local;
2433  if(global.data() == loc.data())
2434  {
2435  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2436  }
2437  else
2438  {
2439  local = loc; // create reference
2440  }
2441 
2442 
2443  if(m_signChange)
2444  {
2445  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2446  }
2447  else
2448  {
2449  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2450  }
2451 
2452  // ensure all values are unique by calling a max
2453  Gs::Gather(global, Gs::gs_max, m_gsh);
2454  }
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
Gs::Gather
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:224
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:115
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:659
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:117
void Nektar::MultiRegions::AssemblyMapCG::v_LocalToGlobal ( const NekVector< NekDouble > &  loc,
NekVector< NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2456 of file AssemblyMapCG.cpp.

References Nektar::NekVector< DataType >::GetPtr(), and Nektar::MultiRegions::AssemblyMap::LocalToGlobal().

2459  {
2460  LocalToGlobal(loc.GetPtr(),global.GetPtr());
2461  }
Nektar::MultiRegions::AssemblyMap::LocalToGlobal
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Definition: AssemblyMap.cpp:722
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2530 of file AssemblyMapCG.cpp.

References Gs::Gather(), Gs::gs_add, and Nektar::MultiRegions::AssemblyMap::m_gsh.

2532  {
2533  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2534  }
Gs::Gather
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:224
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( NekVector< NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2536 of file AssemblyMapCG.cpp.

References Nektar::NekVector< DataType >::GetPtr(), and Nektar::MultiRegions::AssemblyMap::UniversalAssemble().

2538  {
2539  UniversalAssemble(pGlobal.GetPtr());
2540  }
Nektar::MultiRegions::AssemblyMap::UniversalAssemble
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:764
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal,
int  offset 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2542 of file AssemblyMapCG.cpp.

References Nektar::MultiRegions::AssemblyMap::UniversalAssemble(), and Vmath::Vcopy().

2545  {
2546  Array<OneD, NekDouble> tmp(offset);
2547  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2548  UniversalAssemble(pGlobal);
2549  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2550  }
Nektar::MultiRegions::AssemblyMap::UniversalAssemble
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:764
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047

Member Data Documentation

std::map<int, std::vector<ExtraDirDof> > Nektar::MultiRegions::AssemblyMapCG::m_extraDirDofs
protected

Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor.

Definition at line 148 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), and GetExtraDirDofs().

Array<OneD, int> Nektar::MultiRegions::AssemblyMapCG::m_extraDirEdges
protected

Extra dirichlet edges in parallel.

Definition at line 141 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetExtraDirEdges().

int Nektar::MultiRegions::AssemblyMapCG::m_fullSystemBandWidth
protected

Bandwith of the full matrix system (no static condensation).

Definition at line 119 of file AssemblyMapCG.h.

Referenced by CalculateFullSystemBandWidth(), and v_GetFullSystemBandWidth().

Array<OneD,int> Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
protected

Integer map of process coeffs to universal space.

Definition at line 121 of file AssemblyMapCG.h.

Referenced by Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), SetUpUniversalC0ContMap(), v_GetGlobalToUniversalMap(), and v_LinearSpaceMap().

Array<OneD,int> Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMapUnique
protected

Integer map of unique process coeffs to universal space (signed)

Definition at line 123 of file AssemblyMapCG.h.

Referenced by Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), SetUpUniversalC0ContMap(), and v_GetGlobalToUniversalMapUnique().

Array<OneD,int> Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
protected

Integer map of local coeffs to global space.

Definition at line 115 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), CalculateFullSystemBandWidth(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(), SetUpUniversalC0ContMap(), v_Assemble(), v_GetLocalToGlobalMap(), v_GlobalToLocal(), v_LinearSpaceMap(), and v_LocalToGlobal().

Array<OneD,NekDouble> Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
protected

Integer sign of local coeffs to global space.

Definition at line 117 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(), v_Assemble(), v_GetLocalToGlobalSign(), v_GlobalToLocal(), and v_LocalToGlobal().

int Nektar::MultiRegions::AssemblyMapCG::m_maxStaticCondLevel
protected

Maximum static condensation level.

Definition at line 145 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG().

int Nektar::MultiRegions::AssemblyMapCG::m_numDirEdges
protected

Number of Dirichlet edges.

Definition at line 131 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirEdges().

int Nektar::MultiRegions::AssemblyMapCG::m_numDirFaces
protected

Number of Dirichlet faces.

Definition at line 133 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirFaces().

int Nektar::MultiRegions::AssemblyMapCG::m_numLocalBndCondCoeffs
protected

Number of local boundary condition coefficients.

Definition at line 139 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and CreateGraph().

int Nektar::MultiRegions::AssemblyMapCG::m_numLocDirBndCondDofs
protected

Number of local boundary condition degrees of freedom.

Definition at line 143 of file AssemblyMapCG.h.

int Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdgeModes
protected

Number of non Dirichlet edge modes.

Definition at line 127 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().

int Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdges
protected

Number of Dirichlet edges.

Definition at line 135 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdges().

int Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaceModes
protected

Number of non Dirichlet face modes.

Definition at line 129 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().

int Nektar::MultiRegions::AssemblyMapCG::m_numNonDirFaces
protected

Number of Dirichlet faces.

Definition at line 137 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaces().

int Nektar::MultiRegions::AssemblyMapCG::m_numNonDirVertexModes
protected

Number of non Dirichlet vertex modes.

Definition at line 125 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().

Generated on Sat Jun 4 2016 00:02:39 for Nektar++ by   doxygen 1.8.8