Nektar++
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...
 
map< int, vector< ExtraDirDof > > & GetExtraDirDofs ()
 
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
 AssemblyMap ()
 Default constructor. More...
 
 AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Constructor with a communicator. More...
 
 AssemblyMap (AssemblyMap *oldLevelMap, const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
 Constructor for next level in multi-level static condensation. More...
 
virtual ~AssemblyMap ()
 Destructor. More...
 
LibUtilities::CommSharedPtr GetComm ()
 Retrieves the communicator. More...
 
size_t GetHash () const
 Retrieves the hash of this map. More...
 
int GetLocalToGlobalMap (const int i) const
 
int GetGlobalToUniversalMap (const int i) const
 
int GetGlobalToUniversalMapUnique (const int i) const
 
const Array< OneD, const int > & GetLocalToGlobalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMapUnique ()
 
NekDouble GetLocalToGlobalSign (const int i) const
 
const Array< OneD, NekDouble > & GetLocalToGlobalSign () const
 
void LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
void GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
void GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
 
void Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
 
void UniversalAssemble (NekVector< NekDouble > &pGlobal) const
 
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
 
int GetLocalToGlobalBndMap (const int i) const
 Retrieve the global index of a given local boundary mode. More...
 
const Array< OneD, const int > & GetLocalToGlobalBndMap ()
 Retrieve the global indices of the local boundary modes. More...
 
const Array< OneD, const int > & GetGlobalToUniversalBndMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalBndMapUnique ()
 
bool GetSignChange ()
 Returns true if using a modal expansion requiring a change of sign of some modes. More...
 
NekDouble GetLocalToGlobalBndSign (const int i) const
 Retrieve the sign change of a given local boundary mode. More...
 
Array< OneD, const NekDoubleGetLocalToGlobalBndSign () 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< AssemblyMapLinearSpaceMap (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 PatchMapSharedPtrGetPatchMapFromPrevLevel (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 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, set< int > &extraDirVerts, set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
 
void SetUpUniversalC0ContMap (const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
 
void CalculateFullSystemBandWidth ()
 Calculate the bandwith of the full matrix system. More...
 
virtual int v_GetLocalToGlobalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMapUnique (const int i) const
 
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap ()
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap ()
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique ()
 
virtual NekDouble v_GetLocalToGlobalSign (const int i) const
 
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign () const
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
virtual void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
 
virtual void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
virtual void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
 
virtual void v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
 
virtual int v_GetFullSystemBandWidth () const
 
virtual int v_GetNumNonDirVertexModes () const
 
virtual int v_GetNumNonDirEdgeModes () const
 
virtual int v_GetNumNonDirFaceModes () const
 
virtual int v_GetNumDirEdges () const
 
virtual int v_GetNumDirFaces () const
 
virtual int v_GetNumNonDirEdges () const
 
virtual int v_GetNumNonDirFaces () const
 
virtual const Array< OneD, const int > & v_GetExtraDirEdges ()
 
virtual AssemblyMapSharedPtr v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
 Construct an AssemblyMapCG object which corresponds to the linear space of the current object. More...
 
- Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMap
void CalculateBndSystemBandWidth ()
 Calculates the bandwidth of the boundary system. More...
 
void GlobalToLocalBndWithoutSign (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc)
 

Protected Attributes

Array< OneD, int > m_localToGlobalMap
 Integer map of local coeffs to global space. More...
 
Array< OneD, NekDoublem_localToGlobalSign
 Integer sign of local coeffs to global space. More...
 
int m_fullSystemBandWidth
 Bandwith of the full matrix system (no static condensation). More...
 
Array< OneD, int > m_globalToUniversalMap
 Integer map of process coeffs to universal space. More...
 
Array< OneD, int > m_globalToUniversalMapUnique
 Integer map of unique process coeffs to universal space (signed) More...
 
int m_numNonDirVertexModes
 Number of non Dirichlet vertex modes. More...
 
int m_numNonDirEdgeModes
 Number of non Dirichlet edge modes. More...
 
int m_numNonDirFaceModes
 Number of non Dirichlet face modes. More...
 
int m_numDirEdges
 Number of Dirichlet edges. More...
 
int m_numDirFaces
 Number of Dirichlet faces. More...
 
int m_numNonDirEdges
 Number of Dirichlet edges. More...
 
int m_numNonDirFaces
 Number of Dirichlet faces. More...
 
int m_numLocalBndCondCoeffs
 Number of local boundary condition coefficients. More...
 
Array< OneD, int > m_extraDirEdges
 Extra dirichlet edges in parallel. More...
 
int m_numLocDirBndCondDofs
 Number of local boundary condition degrees of freedom. More...
 
int m_maxStaticCondLevel
 Maximum static condensation level. More...
 
map< int, vector< ExtraDirDof > > m_extraDirDofs
 Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. More...
 
- Protected Attributes inherited from Nektar::MultiRegions::AssemblyMap
LibUtilities::SessionReaderSharedPtr m_session
 Session object. More...
 
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
size_t m_hash
 Hash for map. More...
 
int m_numLocalBndCoeffs
 Number of local boundary coefficients. More...
 
int m_numGlobalBndCoeffs
 Total number of global boundary coefficients. More...
 
int m_numLocalDirBndCoeffs
 Number of Local Dirichlet Boundary Coefficients. More...
 
int m_numGlobalDirBndCoeffs
 Number of Global Dirichlet Boundary Coefficients. More...
 
bool m_systemSingular
 Flag indicating if the system is singular or not. More...
 
int m_numLocalCoeffs
 Total number of local coefficients. More...
 
int m_numGlobalCoeffs
 Total number of global coefficients. More...
 
bool m_signChange
 Flag indicating if modes require sign reversal. More...
 
Array< OneD, int > m_localToGlobalBndMap
 Integer map of local boundary coeffs to global space. More...
 
Array< OneD, NekDoublem_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, NekDoublem_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...
 
NekDouble m_iterativeTolerance
 Tolerance for iterative solver. More...
 
int m_successiveRHS
 sucessive RHS for iterative solver More...
 
Gs::gs_datam_gsh
 
Gs::gs_datam_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 ExpListSharedPtrBndCondExp
 
typedef Array< OneD, const SpatialDomains::BoundaryConditionShPtrBndCond
 

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

Definition at line 75 of file AssemblyMapCG.h.

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 69 of file AssemblyMapCG.cpp.

References m_maxStaticCondLevel.

71  :
72  AssemblyMap(pSession,variable)
73  {
74  pSession->LoadParameter(
75  "MaxStaticCondLevel", m_maxStaticCondLevel, 100);
76  }
int m_maxStaticCondLevel
Maximum static condensation level.
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:77
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 1237 of file AssemblyMapCG.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), 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, Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), 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().

1248  : AssemblyMap(pSession, variable)
1249  {
1250  int i, j, k, l;
1251  int cnt = 0;
1252  int intDofCnt;
1253  int meshVertId, meshEdgeId, meshFaceId;
1254  int globalId;
1255  int nEdgeInteriorCoeffs;
1256  int nFaceInteriorCoeffs;
1257  int firstNonDirGraphVertId;
1258  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1261  StdRegions::Orientation edgeOrient;
1262  StdRegions::Orientation faceOrient;
1263  Array<OneD, unsigned int> edgeInteriorMap;
1264  Array<OneD, int> edgeInteriorSign;
1265  Array<OneD, unsigned int> faceInteriorMap;
1266  Array<OneD, int> faceInteriorSign;
1267  PeriodicMap::const_iterator pIt;
1268 
1269  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1270 
1271  m_signChange = false;
1272 
1273  // Stores vertex, edge and face reordered vertices.
1274  DofGraph graph(3);
1275  DofGraph dofs(3);
1276 
1277  set<int> extraDirVerts, extraDirEdges;
1279 
1280  // Construct list of number of degrees of freedom for each vertex,
1281  // edge and face.
1282  for (i = 0; i < locExpVector.size(); ++i)
1283  {
1284  exp = locExpVector[i];
1285 
1286  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
1287  {
1288  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1289  }
1290 
1291  for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
1292  {
1293  dofs[1][exp->GetGeom()->GetEid(j)] =
1294  exp->GetEdgeNcoeffs(j) - 2;
1295  }
1296 
1297  for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
1298  {
1299  dofs[2][exp->GetGeom()->GetFid(j)] =
1300  exp->GetFaceIntNcoeffs(j);
1301  }
1302  }
1303 
1304  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1305 
1306  // Note that nExtraDirichlet is not used in the logic below; it just
1307  // needs to be set so that the coupled solver in
1308  // IncNavierStokesSolver can work.
1309  int nExtraDirichlet;
1310  int nGraphVerts =
1311  CreateGraph(locExp, bndCondExp, bndCondVec,
1312  checkIfSystemSingular, periodicVerts, periodicEdges,
1313  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1314  extraDirEdges, firstNonDirGraphVertId,
1315  nExtraDirichlet);
1316 
1317  /*
1318  * Set up an array which contains the offset information of the
1319  * different graph vertices.
1320  *
1321  * This basically means to identify to how many global degrees of
1322  * freedom the individual graph vertices correspond. Obviously,
1323  * the graph vertices corresponding to the mesh-vertices account
1324  * for a single global DOF. However, the graph vertices
1325  * corresponding to the element edges correspond to N-2 global DOF
1326  * where N is equal to the number of boundary modes on this edge.
1327  */
1328  Array<OneD, int> graphVertOffset(
1329  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1330 
1331  graphVertOffset[0] = 0;
1332 
1333  for(i = 0; i < locExpVector.size(); ++i)
1334  {
1335  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1336 
1337  for(j = 0; j < exp->GetNverts(); ++j)
1338  {
1339  meshVertId = exp->GetGeom()->GetVid(j);
1340  graphVertOffset[graph[0][meshVertId]+1] = 1;
1341  }
1342 
1343  for(j = 0; j < exp->GetNedges(); ++j)
1344  {
1345  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1346  meshEdgeId = exp->GetGeom()->GetEid(j);
1347  graphVertOffset[graph[1][meshEdgeId]+1]
1348  = nEdgeInteriorCoeffs;
1349 
1350  bType = exp->GetEdgeBasisType(j);
1351 
1352  // need a sign vector for modal expansions if nEdgeCoeffs
1353  // >=4
1354  if(nEdgeInteriorCoeffs >= 2 &&
1355  (bType == LibUtilities::eModified_A ||
1356  bType == LibUtilities::eModified_B))
1357  {
1358  m_signChange = true;
1359  }
1360  }
1361 
1362  for(j = 0; j < exp->GetNfaces(); ++j)
1363  {
1364  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1365  meshFaceId = exp->GetGeom()->GetFid(j);
1366  graphVertOffset[graph[2][meshFaceId]+1] = nFaceInteriorCoeffs;
1367  }
1368  }
1369 
1370  for(i = 1; i < graphVertOffset.num_elements(); i++)
1371  {
1372  graphVertOffset[i] += graphVertOffset[i-1];
1373  }
1374 
1375  // Allocate the proper amount of space for the class-data
1376  m_numLocalCoeffs = numLocalCoeffs;
1377  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1381  // If required, set up the sign-vector
1382  if(m_signChange)
1383  {
1387  }
1388 
1389  m_staticCondLevel = 0;
1390  m_numPatches = locExpVector.size();
1393  for(i = 0; i < m_numPatches; ++i)
1394  {
1395  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1396  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1397  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1398  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1399  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1400  }
1401 
1402  /**
1403  * STEP 6: Now, all ingredients are ready to set up the actual
1404  * local to global mapping.
1405  *
1406  * The remainder of the map consists of the element-interior
1407  * degrees of freedom. This leads to the block-diagonal submatrix
1408  * as each element-interior mode is globally orthogonal to modes
1409  * in all other elements.
1410  */
1411  cnt = 0;
1412 
1413  // Loop over all the elements in the domain
1414  for(i = 0; i < locExpVector.size(); ++i)
1415  {
1416  exp = locExpVector[i];
1417  cnt = locExp.GetCoeff_Offset(i);
1418  for(j = 0; j < exp->GetNverts(); ++j)
1419  {
1420  meshVertId = exp->GetGeom()->GetVid(j);
1421 
1422  // Set the global DOF for vertex j of element i
1423  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1424  graphVertOffset[graph[0][meshVertId]];
1425  }
1426 
1427  for(j = 0; j < exp->GetNedges(); ++j)
1428  {
1429  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1430  edgeOrient = exp->GetGeom()->GetEorient(j);
1431  meshEdgeId = exp->GetGeom()->GetEid(j);
1432 
1433  pIt = periodicEdges.find(meshEdgeId);
1434 
1435  // See if this edge is periodic. If it is, then we map all
1436  // edges to the one with lowest ID, and align all
1437  // coefficients to this edge orientation.
1438  if (pIt != periodicEdges.end())
1439  {
1440  pair<int, StdRegions::Orientation> idOrient =
1442  meshEdgeId, edgeOrient, pIt->second);
1443  edgeOrient = idOrient.second;
1444  }
1445 
1446  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1447 
1448  // Set the global DOF's for the interior modes of edge j
1449  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1450  {
1451  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1452  graphVertOffset[graph[1][meshEdgeId]]+k;
1453  }
1454 
1455  // Fill the sign vector if required
1456  if(m_signChange)
1457  {
1458  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1459  {
1460  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1461  }
1462  }
1463  }
1464 
1465  for(j = 0; j < exp->GetNfaces(); ++j)
1466  {
1467  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1468  faceOrient = exp->GetGeom()->GetForient(j);
1469  meshFaceId = exp->GetGeom()->GetFid(j);
1470 
1471  pIt = periodicFaces.find(meshFaceId);
1472 
1473  if (pIt != periodicFaces.end() &&
1474  meshFaceId == min(meshFaceId, pIt->second[0].id))
1475  {
1476  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1477  }
1478 
1479  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1480 
1481  // Set the global DOF's for the interior modes of face j
1482  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1483  {
1484  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1485  graphVertOffset[graph[2][meshFaceId]]+k;
1486  }
1487 
1488  if(m_signChange)
1489  {
1490  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1491  {
1492  m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
1493  }
1494  }
1495  }
1496  }
1497 
1498  // Set up the mapping for the boundary conditions
1499  cnt = 0;
1500  int offset = 0;
1501  for(i = 0; i < bndCondExp.num_elements(); i++)
1502  {
1503  set<int> foundExtraVerts, foundExtraEdges;
1504  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1505  {
1506  bndExp = bndCondExp[i]->GetExp(j);
1507  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1508  for(k = 0; k < bndExp->GetNverts(); k++)
1509  {
1510  meshVertId = bndExp->GetGeom()->GetVid(k);
1511  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1512 
1513  if (bndConditions[i]->GetBoundaryConditionType() !=
1515  {
1516  continue;
1517  }
1518 
1519  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1520  if (iter != extraDirVerts.end() &&
1521  foundExtraVerts.count(meshVertId) == 0)
1522  {
1523  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1524  bndExp->GetVertexMap(k);
1525  int gid = graphVertOffset[
1526  graph[0][meshVertId]];
1527  ExtraDirDof t(loc, gid, 1.0);
1528  m_extraDirDofs[i].push_back(t);
1529  foundExtraVerts.insert(meshVertId);
1530  }
1531  }
1532 
1533  for(k = 0; k < bndExp->GetNedges(); k++)
1534  {
1535  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1536  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1537  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1538 
1539  pIt = periodicEdges.find(meshEdgeId);
1540 
1541  // See if this edge is periodic. If it is, then we map
1542  // all edges to the one with lowest ID, and align all
1543  // coefficients to this edge orientation.
1544  if (pIt != periodicEdges.end())
1545  {
1546  pair<int, StdRegions::Orientation> idOrient =
1548  meshEdgeId, edgeOrient, pIt->second);
1549  edgeOrient = idOrient.second;
1550  }
1551 
1552  bndExp->GetEdgeInteriorMap(
1553  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1554 
1555  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1556  {
1557  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1558  graphVertOffset[graph[1][meshEdgeId]]+l;
1559  }
1560 
1561  // Fill the sign vector if required
1562  if(m_signChange)
1563  {
1564  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1565  {
1566  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1567  }
1568  }
1569 
1570  if (bndConditions[i]->GetBoundaryConditionType() !=
1572  {
1573  continue;
1574  }
1575 
1576  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1577  if (iter != extraDirEdges.end() &&
1578  foundExtraEdges.count(meshEdgeId) == 0 &&
1579  nEdgeInteriorCoeffs > 0)
1580  {
1581  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1582  {
1583  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1584  edgeInteriorMap[l];
1585  int gid = graphVertOffset[
1586  graph[1][meshEdgeId]]+l;
1587  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1588  m_extraDirDofs[i].push_back(t);
1589  }
1590  foundExtraEdges.insert(meshEdgeId);
1591  }
1592  }
1593 
1594  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1595  intDofCnt = 0;
1596  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1597  {
1598  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1599  {
1601  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1602  intDofCnt++;
1603  }
1604  }
1605  }
1606  offset += bndCondExp[i]->GetNcoeffs();
1607  }
1608 
1609  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1610  m_numGlobalBndCoeffs = globalId;
1611 
1612 
1613  /*
1614  * The boundary condition mapping is generated from the same vertex
1615  * renumbering.
1616  */
1617  cnt=0;
1618  for(i = 0; i < m_numLocalCoeffs; ++i)
1619  {
1620  if(m_localToGlobalMap[i] == -1)
1621  {
1622  m_localToGlobalMap[i] = globalId++;
1623  }
1624  else
1625  {
1626  if(m_signChange)
1627  {
1629  }
1631  }
1632  }
1633 
1634  m_numGlobalCoeffs = globalId;
1635 
1636  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1637 
1638  // Since we can now have multiple entries to m_extraDirDofs due to
1639  // periodic boundary conditions we make a call to work out the
1640  // multiplicity of all entries and invert (Need to be after
1641  // SetupUniversalC0ContMap)
1643 
1644  // Fill in Dirichlet coefficients that are to be sent to other
1645  // processors with a value of 1
1646  map<int, vector<ExtraDirDof> >::iterator Tit;
1647 
1648  // Generate valence for extraDirDofs
1649  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1650  {
1651  for (i = 0; i < Tit->second.size(); ++i)
1652  {
1653  valence[Tit->second[i].get<1>()] = 1.0;
1654  }
1655  }
1656 
1657  UniversalAssembleBnd(valence);
1658 
1659  // Set third argument of tuple to inverse of valence.
1660  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1661  {
1662  for (i = 0; i < Tit->second.size(); ++i)
1663  {
1664  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1665  }
1666  }
1667 
1668  // Set up the local to global map for the next level when using
1669  // multi-level static condensation
1673  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1674  {
1675  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1676  {
1677  Array<OneD, int> vwgts_perm(
1678  dofs[0].size() + dofs[1].size() + dofs[2].size()
1679  - firstNonDirGraphVertId);
1680 
1681  for (i = 0; i < locExpVector.size(); ++i)
1682  {
1683  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1684 
1685  for (j = 0; j < exp->GetNverts(); ++j)
1686  {
1687  meshVertId = exp->GetGeom()->GetVid(j);
1688 
1689  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1690  {
1691  vwgts_perm[graph[0][meshVertId] -
1692  firstNonDirGraphVertId] =
1693  dofs[0][meshVertId];
1694  }
1695  }
1696 
1697  for (j = 0; j < exp->GetNedges(); ++j)
1698  {
1699  meshEdgeId = exp->GetGeom()->GetEid(j);
1700 
1701  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
1702  {
1703  vwgts_perm[graph[1][meshEdgeId] -
1704  firstNonDirGraphVertId] =
1705  dofs[1][meshEdgeId];
1706  }
1707  }
1708 
1709  for (j = 0; j < exp->GetNfaces(); ++j)
1710  {
1711  meshFaceId = exp->GetGeom()->GetFid(j);
1712 
1713  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
1714  {
1715  vwgts_perm[graph[2][meshFaceId] -
1716  firstNonDirGraphVertId] =
1717  dofs[2][meshFaceId];
1718  }
1719  }
1720  }
1721 
1722  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1724  AllocateSharedPtr(this, bottomUpGraph);
1725  }
1726  }
1727 
1728  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
1729  m_localToGlobalMap.end());
1730 
1731  // Add up hash values if parallel
1732  int hash = m_hash;
1733  vComm->AllReduce(hash, LibUtilities::ReduceSum);
1734  m_hash = hash;
1735 
1738  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:344
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:313
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:305
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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:756
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
Principle Modified Functions .
Definition: BasisType.h:49
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, set< int > &extraDirVerts, set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:330
vector< map< int, int > > DofGraph
Definition: AssemblyMapCG.h:56
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
void CalculateFullSystemBandWidth()
Calculate the bandwith of the full matrix system.
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:391
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:308
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
map< int, vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:384
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:362
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:353
boost::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:54
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:317
Principle Modified Functions .
Definition: BasisType.h:50
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
double NekDouble
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...
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:386
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
Definition: AssemblyMap.h:347
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:351
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:311
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:380
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:349
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:77
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void UniversalAssembleBnd(Array< OneD, NekDouble > &pGlobal) const
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:341
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:382
Nektar::MultiRegions::AssemblyMapCG::~AssemblyMapCG ( )
virtual

Destructor.

Definition at line 1743 of file AssemblyMapCG.cpp.

1744  {
1745  }

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 2201 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().

2202  {
2203  int i,j;
2204  int cnt = 0;
2205  int locSize;
2206  int maxId;
2207  int minId;
2208  int bwidth = -1;
2209  for(i = 0; i < m_numPatches; ++i)
2210  {
2212  maxId = -1;
2213  minId = m_numLocalCoeffs+1;
2214  for(j = 0; j < locSize; j++)
2215  {
2217  {
2218  if(m_localToGlobalMap[cnt+j] > maxId)
2219  {
2220  maxId = m_localToGlobalMap[cnt+j];
2221  }
2222 
2223  if(m_localToGlobalMap[cnt+j] < minId)
2224  {
2225  minId = m_localToGlobalMap[cnt+j];
2226  }
2227  }
2228  }
2229  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2230 
2231  cnt+=locSize;
2232  }
2233 
2234  m_fullSystemBandWidth = bwidth;
2235  }
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:330
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:384
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:317
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:386
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:382
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,
set< int > &  extraDirVerts,
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 78 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::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().

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

Definition at line 108 of file AssemblyMapCG.h.

References m_extraDirDofs.

109  {
110  return m_extraDirDofs;
111  }
map< int, vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
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 1857 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().

1862  {
1864  int nVert = 0;
1865  int nEdge = 0;
1866  int nFace = 0;
1867  int maxEdgeDof = 0;
1868  int maxFaceDof = 0;
1869  int maxIntDof = 0;
1870  int dof = 0;
1871  int cnt;
1872  int i,j,k;
1873  int meshVertId;
1874  int meshEdgeId;
1875  int meshFaceId;
1876  int elementId;
1877  int vGlobalId;
1878  int maxBndGlobalId = 0;
1879  StdRegions::Orientation edgeOrient;
1880  StdRegions::Orientation faceOrient;
1881  Array<OneD, unsigned int> edgeInteriorMap;
1882  Array<OneD, int> edgeInteriorSign;
1883  Array<OneD, unsigned int> faceInteriorMap;
1884  Array<OneD, int> faceInteriorSign;
1885  Array<OneD, unsigned int> interiorMap;
1886  PeriodicMap::const_iterator pIt;
1887 
1888  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1889  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
1890 
1895 
1896  // Loop over all the elements in the domain to gather mesh data
1897  for(i = 0; i < locExpVector.size(); ++i)
1898  {
1899  exp = locExpVector[i];
1900  nVert += exp->GetNverts();
1901  nEdge += exp->GetNedges();
1902  nFace += exp->GetNfaces();
1903  // Loop over all edges (and vertices) of element i
1904  for(j = 0; j < exp->GetNedges(); ++j)
1905  {
1906  dof = exp->GetEdgeNcoeffs(j)-2;
1907  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
1908  }
1909  for(j = 0; j < exp->GetNfaces(); ++j)
1910  {
1911  dof = exp->GetFaceIntNcoeffs(j);
1912  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
1913  }
1914  exp->GetInteriorMap(interiorMap);
1915  dof = interiorMap.num_elements();
1916  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
1917  }
1918 
1919  // Tell other processes about how many dof we have
1920  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
1921  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
1922  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
1923  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
1924  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
1925  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
1926 
1927  // Assemble global to universal mapping for this process
1928  for(i = 0; i < locExpVector.size(); ++i)
1929  {
1930  exp = locExpVector[i];
1931  cnt = locExp.GetCoeff_Offset(i);
1932 
1933  // Loop over all vertices of element i
1934  for(j = 0; j < exp->GetNverts(); ++j)
1935  {
1936  meshVertId = exp->GetGeom()->GetVid(j);
1937  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
1938 
1939  pIt = perVerts.find(meshVertId);
1940  if (pIt != perVerts.end())
1941  {
1942  for (k = 0; k < pIt->second.size(); ++k)
1943  {
1944  meshVertId = min(meshVertId, pIt->second[k].id);
1945  }
1946  }
1947 
1948  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
1949  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
1950  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
1951  }
1952 
1953  // Loop over all edges of element i
1954  for(j = 0; j < exp->GetNedges(); ++j)
1955  {
1956  meshEdgeId = exp->GetGeom()->GetEid(j);
1957  pIt = perEdges.find(meshEdgeId);
1958  edgeOrient = exp->GetGeom()->GetEorient(j);
1959 
1960  if (pIt != perEdges.end())
1961  {
1962  pair<int, StdRegions::Orientation> idOrient =
1964  meshEdgeId, edgeOrient, pIt->second);
1965  meshEdgeId = idOrient.first;
1966  edgeOrient = idOrient.second;
1967  }
1968 
1969  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1970  dof = exp->GetEdgeNcoeffs(j)-2;
1971 
1972  // Set the global DOF's for the interior modes of edge j
1973  for(k = 0; k < dof; ++k)
1974  {
1975  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
1976  m_globalToUniversalMap[vGlobalId]
1977  = nVert + meshEdgeId * maxEdgeDof + k + 1;
1978  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
1979  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
1980  }
1981  }
1982 
1983  // Loop over all faces of element i
1984  for(j = 0; j < exp->GetNfaces(); ++j)
1985  {
1986  faceOrient = exp->GetGeom()->GetForient(j);
1987 
1988  meshFaceId = exp->GetGeom()->GetFid(j);
1989 
1990  pIt = perFaces.find(meshFaceId);
1991  if (pIt != perFaces.end())
1992  {
1993  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
1994  {
1995  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1996  }
1997  meshFaceId = min(meshFaceId, pIt->second[0].id);
1998  }
1999 
2000 
2001  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2002  dof = exp->GetFaceIntNcoeffs(j);
2003 
2004 
2005  for(k = 0; k < dof; ++k)
2006  {
2007  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2008  m_globalToUniversalMap[vGlobalId]
2009  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2010  + k + 1;
2011  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2012 
2013  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2014  }
2015  }
2016 
2017  // Add interior DOFs to complete universal numbering
2018  exp->GetInteriorMap(interiorMap);
2019  dof = interiorMap.num_elements();
2020  elementId = (exp->GetGeom())->GetGlobalID();
2021  for (k = 0; k < dof; ++k)
2022  {
2023  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2024  m_globalToUniversalMap[vGlobalId]
2025  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2026  }
2027  }
2028 
2029  // Set up the GSLib universal assemble mapping
2030  // Internal DOF do not participate in any data
2031  // exchange, so we keep these set to the special GSLib id=0 so
2032  // they are ignored.
2033  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2034  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2035  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2036  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2037  {
2038  tmp[i] = m_globalToUniversalMap[i];
2039  }
2040 
2041  m_gsh = Gs::Init(tmp, vCommRow);
2042  m_bndGsh = Gs::Init(tmp2, vCommRow);
2043  Gs::Unique(tmp, vCommRow);
2044  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2045  {
2046  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2047  }
2048  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2049  {
2050  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2051  }
2052  }
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:313
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:305
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:150
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
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:180
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...
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:357
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:359
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:341
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 2356 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().

2359  {
2361  if(global.data() == loc.data())
2362  {
2363  local = Array<OneD, NekDouble>(local.num_elements(),local.data());
2364  }
2365  else
2366  {
2367  local = loc; // create reference
2368  }
2369  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2370 
2371  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2372 
2373  if(m_signChange)
2374  {
2375  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2376  }
2377  else
2378  {
2379  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2380  }
2381  UniversalAssemble(global);
2382  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:344
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:330
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
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:686
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:341
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 2384 of file AssemblyMapCG.cpp.

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

2387  {
2388  Assemble(loc.GetPtr(),global.GetPtr());
2389  }
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetExtraDirEdges ( )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2453 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2454  {
2455  return m_extraDirEdges;
2456  }
Array< OneD, int > m_extraDirEdges
Extra dirichlet edges in parallel.
int Nektar::MultiRegions::AssemblyMapCG::v_GetFullSystemBandWidth ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2413 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2414  {
2415  return m_fullSystemBandWidth;
2416  }
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2243 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2244  {
2245  return m_globalToUniversalMap[i];
2246  }
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2260 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2261  {
2262  return m_globalToUniversalMap;
2263  }
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2248 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2249  {
2250  return m_globalToUniversalMapUnique[i];
2251  }
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2266 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2267  {
2269  }
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
int Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2238 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2239  {
2240  return m_localToGlobalMap[i];
2241  }
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2254 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2255  {
2256  return m_localToGlobalMap;
2257  }
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
NekDouble Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2271 of file AssemblyMapCG.cpp.

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

2273  {
2274  if(m_signChange)
2275  {
2276  return m_localToGlobalSign[i];
2277  }
2278  else
2279  {
2280  return 1.0;
2281  }
2282  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:344
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
const Array< OneD, NekDouble > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2284 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2285  {
2286  return m_localToGlobalSign;
2287  }
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2433 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2434  {
2435  return m_numDirEdges;
2436  }
int m_numDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2438 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2439  {
2440  return m_numDirFaces;
2441  }
int m_numDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2423 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2424  {
2425  return m_numNonDirEdgeModes;
2426  }
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2443 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2444  {
2445  return m_numNonDirEdges;
2446  }
int m_numNonDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2428 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2429  {
2430  return m_numNonDirFaceModes;
2431  }
int m_numNonDirFaceModes
Number of non Dirichlet face modes.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2448 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2449  {
2450  return m_numNonDirFaces;
2451  }
int m_numNonDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2418 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2419  {
2420  return m_numNonDirVertexModes;
2421  }
int m_numNonDirVertexModes
Number of non Dirichlet vertex modes.
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 2324 of file AssemblyMapCG.cpp.

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

2327  {
2329  if(global.data() == loc.data())
2330  {
2331  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2332  }
2333  else
2334  {
2335  glo = global; // create reference
2336  }
2337 
2338 
2339  if(m_signChange)
2340  {
2341  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2342  }
2343  else
2344  {
2345  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2346  }
2347  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:344
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:330
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
void Nektar::MultiRegions::AssemblyMapCG::v_GlobalToLocal ( const NekVector< NekDouble > &  global,
NekVector< NekDouble > &  loc 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2349 of file AssemblyMapCG.cpp.

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

2352  {
2353  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2354  }
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
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 2062 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().

2064  {
2065  AssemblyMapCGSharedPtr returnval;
2066 
2067  int i, j;
2068  int nverts = 0;
2069  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2070  = locexp.GetExp();
2071  int nelmts = exp->size();
2072 
2073  // Get Default Map and turn off any searched values.
2074  returnval = MemoryManager<AssemblyMapCG>
2076  returnval->m_solnType = solnType;
2077  returnval->m_preconType = eNull;
2078  returnval->m_maxStaticCondLevel = 0;
2079  returnval->m_signChange = false;
2080  returnval->m_comm = m_comm;
2081 
2082  // Count the number of vertices
2083  for (i = 0; i < nelmts; ++i)
2084  {
2085  nverts += (*exp)[i]->GetNverts();
2086  }
2087 
2088  returnval->m_numLocalCoeffs = nverts;
2089  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2090 
2091  // Store original global ids in this map
2092  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2093 
2094  int cnt = 0;
2095  int cnt1 = 0;
2096  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2097 
2098  // Set up local to global map;
2099  for (i = 0; i < nelmts; ++i)
2100  {
2101  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2102  {
2103  returnval->m_localToGlobalMap[cnt] =
2104  returnval->m_localToGlobalBndMap[cnt] =
2105  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2106  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2107 
2108  // Set up numLocalDirBndCoeffs
2109  if ((returnval->m_localToGlobalMap[cnt]) <
2111  {
2112  returnval->m_numLocalDirBndCoeffs++;
2113  }
2114  cnt++;
2115  }
2116  cnt1 += (*exp)[i]->GetNcoeffs();
2117  }
2118 
2119  cnt = 0;
2120  // Reset global numbering and count number of dofs
2121  for (i = 0; i < m_numGlobalCoeffs; ++i)
2122  {
2123  if (GlobCoeffs[i] != -1)
2124  {
2125  GlobCoeffs[i] = cnt++;
2126  }
2127  }
2128 
2129  // Set up number of globalCoeffs;
2130  returnval->m_numGlobalCoeffs = cnt;
2131 
2132  // Set up number of global Dirichlet boundary coefficients
2133  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2134  {
2135  if (GlobCoeffs[i] != -1)
2136  {
2137  returnval->m_numGlobalDirBndCoeffs++;
2138  }
2139  }
2140 
2141  // Set up global to universal map
2142  if (m_globalToUniversalMap.num_elements())
2143  {
2145  = m_session->GetComm()->GetRowComm();
2146  int nglocoeffs = returnval->m_numGlobalCoeffs;
2147  returnval->m_globalToUniversalMap
2148  = Array<OneD, int> (nglocoeffs);
2149  returnval->m_globalToUniversalMapUnique
2150  = Array<OneD, int> (nglocoeffs);
2151 
2152  // Reset local to global map and setup universal map
2153  for (i = 0; i < nverts; ++i)
2154  {
2155  cnt = returnval->m_localToGlobalMap[i];
2156  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2157 
2158  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2160  }
2161 
2162  Nektar::Array<OneD, long> tmp(nglocoeffs);
2163  Vmath::Zero(nglocoeffs, tmp, 1);
2164  for (unsigned int i = 0; i < nglocoeffs; ++i)
2165  {
2166  tmp[i] = returnval->m_globalToUniversalMap[i];
2167  }
2168  returnval->m_gsh = Gs::Init(tmp, vCommRow);
2169  Gs::Unique(tmp, vCommRow);
2170  for (unsigned int i = 0; i < nglocoeffs; ++i)
2171  {
2172  returnval->m_globalToUniversalMapUnique[i]
2173  = (tmp[i] >= 0 ? 1 : 0);
2174  }
2175  }
2176  else // not sure this option is ever needed.
2177  {
2178  for (i = 0; i < nverts; ++i)
2179  {
2180  cnt = returnval->m_localToGlobalMap[i];
2181  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2182  }
2183  }
2184 
2185  return returnval;
2186  }
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:305
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:150
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:317
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:180
No Solution type specified.
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:302
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:341
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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 2289 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().

2292  {
2294  if(global.data() == loc.data())
2295  {
2296  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2297  }
2298  else
2299  {
2300  local = loc; // create reference
2301  }
2302 
2303 
2304  if(m_signChange)
2305  {
2306  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2307  }
2308  else
2309  {
2310  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2311  }
2312 
2313  // ensure all values are unique by calling a max
2314  Gs::Gather(global, Gs::gs_max, m_gsh);
2315  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:344
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:218
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:330
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:659
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
void Nektar::MultiRegions::AssemblyMapCG::v_LocalToGlobal ( const NekVector< NekDouble > &  loc,
NekVector< NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2317 of file AssemblyMapCG.cpp.

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

2320  {
2321  LocalToGlobal(loc.GetPtr(),global.GetPtr());
2322  }
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2391 of file AssemblyMapCG.cpp.

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

2393  {
2394  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2395  }
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:218
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( NekVector< NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2397 of file AssemblyMapCG.cpp.

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

2399  {
2400  UniversalAssemble(pGlobal.GetPtr());
2401  }
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal,
int  offset 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2403 of file AssemblyMapCG.cpp.

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

2406  {
2407  Array<OneD, NekDouble> tmp(offset);
2408  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2409  UniversalAssemble(pGlobal);
2410  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2411  }
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038

Member Data Documentation

map<int, 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
Array<OneD,NekDouble> Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
protected
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().