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 64 of file AssemblyMapCG.h.

Member Typedef Documentation

Definition at line 68 of file AssemblyMapCG.h.

Definition at line 66 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 1235 of file AssemblyMapCG.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::AssemblyMap::CalculateBndSystemBandWidth(), CalculateFullSystemBandWidth(), CreateGraph(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::StdRegions::eBackwards, Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::StdRegions::eForwards, 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().

1246  : AssemblyMap(pSession, variable)
1247  {
1248  int i, j, k, l;
1249  int cnt = 0;
1250  int intDofCnt;
1251  int meshVertId, meshEdgeId, meshFaceId;
1252  int globalId;
1253  int nEdgeInteriorCoeffs;
1254  int nFaceInteriorCoeffs;
1255  int firstNonDirGraphVertId;
1256  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1259  StdRegions::Orientation edgeOrient;
1260  StdRegions::Orientation faceOrient;
1261  Array<OneD, unsigned int> edgeInteriorMap;
1262  Array<OneD, int> edgeInteriorSign;
1263  Array<OneD, unsigned int> faceInteriorMap;
1264  Array<OneD, int> faceInteriorSign;
1265  PeriodicMap::const_iterator pIt;
1266 
1267  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1268 
1269  m_signChange = false;
1270 
1271  // Stores vertex, edge and face reordered vertices.
1272  DofGraph graph(3);
1273  DofGraph dofs(3);
1274 
1275  set<int> extraDirVerts, extraDirEdges;
1277 
1278  // Construct list of number of degrees of freedom for each vertex,
1279  // edge and face.
1280  for (i = 0; i < locExpVector.size(); ++i)
1281  {
1282  exp = locExpVector[i];
1283 
1284  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
1285  {
1286  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1287  }
1288 
1289  for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
1290  {
1291  dofs[1][exp->GetGeom()->GetEid(j)] =
1292  exp->GetEdgeNcoeffs(j) - 2;
1293  }
1294 
1295  for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
1296  {
1297  dofs[2][exp->GetGeom()->GetFid(j)] =
1298  exp->GetFaceIntNcoeffs(j);
1299  }
1300  }
1301 
1302  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1303 
1304  // Note that nExtraDirichlet is not used in the logic below; it just
1305  // needs to be set so that the coupled solver in
1306  // IncNavierStokesSolver can work.
1307  int nExtraDirichlet;
1308  int nGraphVerts =
1309  CreateGraph(locExp, bndCondExp, bndCondVec,
1310  checkIfSystemSingular, periodicVerts, periodicEdges,
1311  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1312  extraDirEdges, firstNonDirGraphVertId,
1313  nExtraDirichlet);
1314 
1315  /*
1316  * Set up an array which contains the offset information of the
1317  * different graph vertices.
1318  *
1319  * This basically means to identify to how many global degrees of
1320  * freedom the individual graph vertices correspond. Obviously,
1321  * the graph vertices corresponding to the mesh-vertices account
1322  * for a single global DOF. However, the graph vertices
1323  * corresponding to the element edges correspond to N-2 global DOF
1324  * where N is equal to the number of boundary modes on this edge.
1325  */
1326  Array<OneD, int> graphVertOffset(
1327  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1328 
1329  graphVertOffset[0] = 0;
1330 
1331  for(i = 0; i < locExpVector.size(); ++i)
1332  {
1333  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1334 
1335  for(j = 0; j < exp->GetNverts(); ++j)
1336  {
1337  meshVertId = exp->GetGeom()->GetVid(j);
1338  graphVertOffset[graph[0][meshVertId]+1] = 1;
1339  }
1340 
1341  for(j = 0; j < exp->GetNedges(); ++j)
1342  {
1343  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1344  meshEdgeId = exp->GetGeom()->GetEid(j);
1345  graphVertOffset[graph[1][meshEdgeId]+1]
1346  = nEdgeInteriorCoeffs;
1347 
1348  bType = exp->GetEdgeBasisType(j);
1349 
1350  // need a sign vector for modal expansions if nEdgeCoeffs
1351  // >=4
1352  if(nEdgeInteriorCoeffs >= 2 &&
1353  (bType == LibUtilities::eModified_A ||
1354  bType == LibUtilities::eModified_B))
1355  {
1356  m_signChange = true;
1357  }
1358  }
1359 
1360  for(j = 0; j < exp->GetNfaces(); ++j)
1361  {
1362  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1363  meshFaceId = exp->GetGeom()->GetFid(j);
1364  graphVertOffset[graph[2][meshFaceId]+1] = nFaceInteriorCoeffs;
1365  }
1366  }
1367 
1368  for(i = 1; i < graphVertOffset.num_elements(); i++)
1369  {
1370  graphVertOffset[i] += graphVertOffset[i-1];
1371  }
1372 
1373  // Allocate the proper amount of space for the class-data
1374  m_numLocalCoeffs = numLocalCoeffs;
1375  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1379  // If required, set up the sign-vector
1380  if(m_signChange)
1381  {
1385  }
1386 
1387  m_staticCondLevel = 0;
1388  m_numPatches = locExpVector.size();
1391  for(i = 0; i < m_numPatches; ++i)
1392  {
1393  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1394  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1395  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1396  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1397  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1398  }
1399 
1400  /**
1401  * STEP 6: Now, all ingredients are ready to set up the actual
1402  * local to global mapping.
1403  *
1404  * The remainder of the map consists of the element-interior
1405  * degrees of freedom. This leads to the block-diagonal submatrix
1406  * as each element-interior mode is globally orthogonal to modes
1407  * in all other elements.
1408  */
1409  cnt = 0;
1410 
1411  // Loop over all the elements in the domain
1412  for(i = 0; i < locExpVector.size(); ++i)
1413  {
1414  exp = locExpVector[i];
1415  cnt = locExp.GetCoeff_Offset(i);
1416  for(j = 0; j < exp->GetNverts(); ++j)
1417  {
1418  meshVertId = exp->GetGeom()->GetVid(j);
1419 
1420  // Set the global DOF for vertex j of element i
1421  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1422  graphVertOffset[graph[0][meshVertId]];
1423  }
1424 
1425  for(j = 0; j < exp->GetNedges(); ++j)
1426  {
1427  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1428  edgeOrient = exp->GetGeom()->GetEorient(j);
1429  meshEdgeId = exp->GetGeom()->GetEid(j);
1430 
1431  pIt = periodicEdges.find(meshEdgeId);
1432 
1433  // See if this edge is periodic. If it is, then we map all
1434  // edges to the one with lowest ID, and align all
1435  // coefficients to this edge orientation.
1436  if (pIt != periodicEdges.end())
1437  {
1438  int minId = pIt->second[0].id;
1439  int minIdK = 0;
1440  for (k = 1; k < pIt->second.size(); ++k)
1441  {
1442  if (pIt->second[k].id < minId)
1443  {
1444  minId = min(minId, pIt->second[k].id);
1445  minIdK = k;
1446  }
1447  }
1448 
1449  if( meshEdgeId != min(minId, meshEdgeId))
1450  {
1451  if (pIt->second[minIdK].orient == StdRegions::eBackwards)
1452  {
1453  // Swap edge orientation
1454  edgeOrient = (edgeOrient == StdRegions::eForwards) ?
1456  }
1457  }
1458 
1459  }
1460 
1461  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1462 
1463  // Set the global DOF's for the interior modes of edge j
1464  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1465  {
1466  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1467  graphVertOffset[graph[1][meshEdgeId]]+k;
1468  }
1469 
1470  // Fill the sign vector if required
1471  if(m_signChange)
1472  {
1473  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1474  {
1475  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1476  }
1477  }
1478  }
1479 
1480  for(j = 0; j < exp->GetNfaces(); ++j)
1481  {
1482  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1483  faceOrient = exp->GetGeom()->GetForient(j);
1484  meshFaceId = exp->GetGeom()->GetFid(j);
1485 
1486  pIt = periodicFaces.find(meshFaceId);
1487 
1488  if (pIt != periodicFaces.end() &&
1489  meshFaceId == min(meshFaceId, pIt->second[0].id))
1490  {
1491  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1492  }
1493 
1494  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1495 
1496  // Set the global DOF's for the interior modes of face j
1497  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1498  {
1499  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1500  graphVertOffset[graph[2][meshFaceId]]+k;
1501  }
1502 
1503  if(m_signChange)
1504  {
1505  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1506  {
1507  m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
1508  }
1509  }
1510  }
1511  }
1512 
1513  // Set up the mapping for the boundary conditions
1514  cnt = 0;
1515  int offset = 0;
1516  for(i = 0; i < bndCondExp.num_elements(); i++)
1517  {
1518  set<int> foundExtraVerts, foundExtraEdges;
1519  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1520  {
1521  bndExp = bndCondExp[i]->GetExp(j);
1522  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1523  for(k = 0; k < bndExp->GetNverts(); k++)
1524  {
1525  meshVertId = bndExp->GetGeom()->GetVid(k);
1526  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1527 
1528  if (bndConditions[i]->GetBoundaryConditionType() !=
1530  {
1531  continue;
1532  }
1533 
1534  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1535  if (iter != extraDirVerts.end() &&
1536  foundExtraVerts.count(meshVertId) == 0)
1537  {
1538  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1539  bndExp->GetVertexMap(k);
1540  int gid = graphVertOffset[
1541  graph[0][meshVertId]];
1542  ExtraDirDof t(loc, gid, 1.0);
1543  m_extraDirDofs[i].push_back(t);
1544  foundExtraVerts.insert(meshVertId);
1545  }
1546  }
1547 
1548  for(k = 0; k < bndExp->GetNedges(); k++)
1549  {
1550  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1551  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1552  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1553 
1554  pIt = periodicEdges.find(meshEdgeId);
1555 
1556  // See if this edge is periodic. If it is, then we map
1557  // all edges to the one with lowest ID, and align all
1558  // coefficients to this edge orientation.
1559  if (pIt != periodicEdges.end())
1560  {
1561  int minId = pIt->second[0].id;
1562  int minIdL = 0;
1563  for (l = 1; l < pIt->second.size(); ++l)
1564  {
1565  if (pIt->second[l].id < minId)
1566  {
1567  minId = min(minId, pIt->second[l].id);
1568  minIdL = l;
1569  }
1570  }
1571 
1572  if (pIt->second[minIdL].orient == StdRegions::eBackwards &&
1573  meshEdgeId != min(minId, meshEdgeId))
1574  {
1575  edgeOrient = (edgeOrient == StdRegions::eForwards) ? StdRegions::eBackwards : StdRegions::eForwards;
1576  }
1577  }
1578 
1579  bndExp->GetEdgeInteriorMap(
1580  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1581 
1582  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1583  {
1584  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1585  graphVertOffset[graph[1][meshEdgeId]]+l;
1586  }
1587 
1588  // Fill the sign vector if required
1589  if(m_signChange)
1590  {
1591  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1592  {
1593  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1594  }
1595  }
1596 
1597  if (bndConditions[i]->GetBoundaryConditionType() !=
1599  {
1600  continue;
1601  }
1602 
1603  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1604  if (iter != extraDirEdges.end() &&
1605  foundExtraEdges.count(meshEdgeId) == 0 &&
1606  nEdgeInteriorCoeffs > 0)
1607  {
1608  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1609  {
1610  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1611  edgeInteriorMap[l];
1612  int gid = graphVertOffset[
1613  graph[1][meshEdgeId]]+l;
1614  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1615  m_extraDirDofs[i].push_back(t);
1616  }
1617  foundExtraEdges.insert(meshEdgeId);
1618  }
1619  }
1620 
1621  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1622  intDofCnt = 0;
1623  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1624  {
1625  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1626  {
1628  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1629  intDofCnt++;
1630  }
1631  }
1632  }
1633  offset += bndCondExp[i]->GetNcoeffs();
1634  }
1635 
1636  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1637  m_numGlobalBndCoeffs = globalId;
1638 
1639 
1640  /*
1641  * The boundary condition mapping is generated from the same vertex
1642  * renumbering.
1643  */
1644  cnt=0;
1645  for(i = 0; i < m_numLocalCoeffs; ++i)
1646  {
1647  if(m_localToGlobalMap[i] == -1)
1648  {
1649  m_localToGlobalMap[i] = globalId++;
1650  }
1651  else
1652  {
1653  if(m_signChange)
1654  {
1656  }
1658  }
1659  }
1660 
1661  m_numGlobalCoeffs = globalId;
1662 
1663  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1664 
1665  // Since we can now have multiple entries to m_extraDirDofs due to
1666  // periodic boundary conditions we make a call to work out the
1667  // multiplicity of all entries and invert (Need to be after
1668  // SetupUniversalC0ContMap)
1670 
1671  // Fill in Dirichlet coefficients that are to be sent to other
1672  // processors with a value of 1
1673  map<int, vector<ExtraDirDof> >::iterator Tit;
1674 
1675  // Generate valence for extraDirDofs
1676  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1677  {
1678  for (i = 0; i < Tit->second.size(); ++i)
1679  {
1680  valence[Tit->second[i].get<1>()] = 1.0;
1681  }
1682  }
1683 
1684  UniversalAssembleBnd(valence);
1685 
1686  // Set third argument of tuple to inverse of valence.
1687  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1688  {
1689  for (i = 0; i < Tit->second.size(); ++i)
1690  {
1691  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1692  }
1693  }
1694 
1695  // Set up the local to global map for the next level when using
1696  // multi-level static condensation
1700  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1701  {
1702  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1703  {
1704  Array<OneD, int> vwgts_perm(
1705  dofs[0].size() + dofs[1].size() + dofs[2].size()
1706  - firstNonDirGraphVertId);
1707 
1708  for (i = 0; i < locExpVector.size(); ++i)
1709  {
1710  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1711 
1712  for (j = 0; j < exp->GetNverts(); ++j)
1713  {
1714  meshVertId = exp->GetGeom()->GetVid(j);
1715 
1716  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1717  {
1718  vwgts_perm[graph[0][meshVertId] -
1719  firstNonDirGraphVertId] =
1720  dofs[0][meshVertId];
1721  }
1722  }
1723 
1724  for (j = 0; j < exp->GetNedges(); ++j)
1725  {
1726  meshEdgeId = exp->GetGeom()->GetEid(j);
1727 
1728  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
1729  {
1730  vwgts_perm[graph[1][meshEdgeId] -
1731  firstNonDirGraphVertId] =
1732  dofs[1][meshEdgeId];
1733  }
1734  }
1735 
1736  for (j = 0; j < exp->GetNfaces(); ++j)
1737  {
1738  meshFaceId = exp->GetGeom()->GetFid(j);
1739 
1740  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
1741  {
1742  vwgts_perm[graph[2][meshFaceId] -
1743  firstNonDirGraphVertId] =
1744  dofs[2][meshFaceId];
1745  }
1746  }
1747  }
1748 
1749  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1751  AllocateSharedPtr(this, bottomUpGraph);
1752  }
1753  }
1754 
1755  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
1756  m_localToGlobalMap.end());
1757 
1758  // Add up hash values if parallel
1759  int hash = m_hash;
1760  vComm->AllReduce(hash, LibUtilities::ReduceSum);
1761  m_hash = hash;
1762 
1765  }
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
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 1770 of file AssemblyMapCG.cpp.

1771  {
1772  }

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

2168  {
2169  int i,j;
2170  int cnt = 0;
2171  int locSize;
2172  int maxId;
2173  int minId;
2174  int bwidth = -1;
2175  for(i = 0; i < m_numPatches; ++i)
2176  {
2178  maxId = -1;
2179  minId = m_numLocalCoeffs+1;
2180  for(j = 0; j < locSize; j++)
2181  {
2183  {
2184  if(m_localToGlobalMap[cnt+j] > maxId)
2185  {
2186  maxId = m_localToGlobalMap[cnt+j];
2187  }
2188 
2189  if(m_localToGlobalMap[cnt+j] < minId)
2190  {
2191  minId = m_localToGlobalMap[cnt+j];
2192  }
2193  }
2194  }
2195  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2196 
2197  cnt+=locSize;
2198  }
2199 
2200  m_fullSystemBandWidth = bwidth;
2201  }
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 
1038  {
1039  vector<long> procVerts, procEdges, procFaces;
1040  set <int> foundVerts, foundEdges, foundFaces;
1041 
1042  // Loop over element and construct the procVerts and procEdges
1043  // vectors, which store the geometry IDs of mesh vertices and
1044  // edges respectively which are local to this process.
1045  for(i = cnt = 0; i < locExpVector.size(); ++i)
1046  {
1047  int elmtid = locExp.GetOffset_Elmt_Id(i);
1048  exp = locExpVector[elmtid];
1049  for (j = 0; j < exp->GetNverts(); ++j)
1050  {
1051  int vid = exp->GetGeom()->GetVid(j)+1;
1052  if (foundVerts.count(vid) == 0)
1053  {
1054  procVerts.push_back(vid);
1055  foundVerts.insert(vid);
1056  }
1057  }
1058 
1059  for (j = 0; j < exp->GetNedges(); ++j)
1060  {
1061  int eid = exp->GetGeom()->GetEid(j)+1;
1062 
1063  if (foundEdges.count(eid) == 0)
1064  {
1065  procEdges.push_back(eid);
1066  foundEdges.insert(eid);
1067  }
1068  }
1069 
1070  for (j = 0; j < exp->GetNfaces(); ++j)
1071  {
1072  int fid = exp->GetGeom()->GetFid(j)+1;
1073 
1074  if (foundFaces.count(fid) == 0)
1075  {
1076  procFaces.push_back(fid);
1077  foundFaces.insert(fid);
1078  }
1079  }
1080  }
1081 
1082  int unique_verts = foundVerts.size();
1083  int unique_edges = foundEdges.size();
1084  int unique_faces = foundFaces.size();
1085 
1086  // Now construct temporary GS objects. These will be used to
1087  // populate the arrays tmp3 and tmp4 with the multiplicity of
1088  // the vertices and edges respectively to identify those
1089  // vertices and edges which are located on partition boundary.
1090  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1091  Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm);
1092  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1093  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1094  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1095  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1096 
1097  if (unique_edges > 0)
1098  {
1099  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1100  Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
1101  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1102  }
1103 
1104  if (unique_faces > 0)
1105  {
1106  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1107  Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
1108  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1109  }
1110 
1111  // Finally, fill the partVerts set with all non-Dirichlet
1112  // vertices which lie on a partition boundary.
1113  for (i = 0; i < unique_verts; ++i)
1114  {
1115  if (tmp4[i] > 1.0)
1116  {
1117  if (graph[0].count(procVerts[i]-1) == 0)
1118  {
1119  partVerts.insert(tempGraph[0][procVerts[i]-1]);
1120  }
1121  }
1122  }
1123 
1124  for (i = 0; i < unique_edges; ++i)
1125  {
1126  if (tmp5[i] > 1.0)
1127  {
1128  if (graph[1].count(procEdges[i]-1) == 0)
1129  {
1130  partVerts.insert(tempGraph[1][procEdges[i]-1]);
1131  }
1132  }
1133  }
1134 
1135  for (i = 0; i < unique_faces; ++i)
1136  {
1137  if (tmp6[i] > 1.0)
1138  {
1139  if (graph[2].count(procFaces[i]-1) == 0)
1140  {
1141  partVerts.insert(tempGraph[2][procFaces[i]-1]);
1142  }
1143  }
1144  }
1145  }
1146 
1147  int nGraphVerts = tempGraphVertId;
1148  Array<OneD, int> perm(nGraphVerts);
1149  Array<OneD, int> iperm(nGraphVerts);
1150  Array<OneD, int> vwgts(nGraphVerts);
1151  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1152  for(i = 0; i < nGraphVerts; ++i)
1153  {
1154  vwgts[i] = vwgts_map[i];
1155  }
1156 
1157  if(nGraphVerts)
1158  {
1159  switch(m_solnType)
1160  {
1161  case eDirectFullMatrix:
1162  case eIterativeFull:
1163  case eIterativeStaticCond:
1164  case ePETScStaticCond:
1165  case ePETScFullMatrix:
1166  case eXxtFullMatrix:
1167  case eXxtStaticCond:
1168  {
1169  NoReordering(boostGraphObj,perm,iperm);
1170  break;
1171  }
1172 
1173  case eDirectStaticCond:
1174  {
1175  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1176  break;
1177  }
1178 
1183  {
1185  boostGraphObj, perm, iperm, bottomUpGraph,
1186  partVerts, mdswitch);
1187  break;
1188  }
1189  default:
1190  {
1191  ASSERTL0(false,
1192  "Unrecognised solution type: " + std::string(
1194  }
1195  }
1196  }
1197 
1198  // For parallel multi-level static condensation determine the lowest
1199  // static condensation level amongst processors.
1201  {
1202  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1203  vComm->AllReduce(m_lowestStaticCondLevel,
1205  }
1206  else
1207  {
1209  }
1210 
1211  /**
1212  * STEP 4: Fill the #graph[0] and
1213  * #graph[1] with the optimal ordering from boost.
1214  */
1216  for(mapIt = tempGraph[0].begin(); mapIt != tempGraph[0].end(); mapIt++)
1217  {
1218  graph[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1219  }
1220  for(mapIt = tempGraph[1].begin(); mapIt != tempGraph[1].end(); mapIt++)
1221  {
1222  graph[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1223  }
1224  for(mapIt = tempGraph[2].begin(); mapIt != tempGraph[2].end(); mapIt++)
1225  {
1226  graph[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1227  }
1228 
1229  return nGraphVerts;
1230  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
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:165
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 101 of file AssemblyMapCG.h.

References m_extraDirDofs.

102  {
103  return m_extraDirDofs;
104  }
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 1825 of file AssemblyMapCG.cpp.

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

1830  {
1832  int nVert = 0;
1833  int nEdge = 0;
1834  int nFace = 0;
1835  int maxEdgeDof = 0;
1836  int maxFaceDof = 0;
1837  int maxIntDof = 0;
1838  int dof = 0;
1839  int cnt;
1840  int i,j,k;
1841  int meshVertId;
1842  int meshEdgeId;
1843  int meshFaceId;
1844  int elementId;
1845  int vGlobalId;
1846  int maxBndGlobalId = 0;
1847  StdRegions::Orientation edgeOrient;
1848  StdRegions::Orientation faceOrient;
1849  Array<OneD, unsigned int> edgeInteriorMap;
1850  Array<OneD, int> edgeInteriorSign;
1851  Array<OneD, unsigned int> faceInteriorMap;
1852  Array<OneD, int> faceInteriorSign;
1853  Array<OneD, unsigned int> interiorMap;
1854  PeriodicMap::const_iterator pIt;
1855 
1856  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1857  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
1858 
1863 
1864  // Loop over all the elements in the domain to gather mesh data
1865  for(i = 0; i < locExpVector.size(); ++i)
1866  {
1867  exp = locExpVector[i];
1868  nVert += exp->GetNverts();
1869  nEdge += exp->GetNedges();
1870  nFace += exp->GetNfaces();
1871  // Loop over all edges (and vertices) of element i
1872  for(j = 0; j < exp->GetNedges(); ++j)
1873  {
1874  dof = exp->GetEdgeNcoeffs(j)-2;
1875  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
1876  }
1877  for(j = 0; j < exp->GetNfaces(); ++j)
1878  {
1879  dof = exp->GetFaceIntNcoeffs(j);
1880  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
1881  }
1882  exp->GetInteriorMap(interiorMap);
1883  dof = interiorMap.num_elements();
1884  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
1885  }
1886 
1887  // Tell other processes about how many dof we have
1888  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
1889  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
1890  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
1891  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
1892  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
1893  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
1894 
1895  // Assemble global to universal mapping for this process
1896  for(i = 0; i < locExpVector.size(); ++i)
1897  {
1898  exp = locExpVector[i];
1899  cnt = locExp.GetCoeff_Offset(i);
1900 
1901  // Loop over all vertices of element i
1902  for(j = 0; j < exp->GetNverts(); ++j)
1903  {
1904  meshVertId = exp->GetGeom()->GetVid(j);
1905  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
1906 
1907  pIt = perVerts.find(meshVertId);
1908  if (pIt != perVerts.end())
1909  {
1910  for (k = 0; k < pIt->second.size(); ++k)
1911  {
1912  meshVertId = min(meshVertId, pIt->second[k].id);
1913  }
1914  }
1915 
1916  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
1917  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
1918  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
1919  }
1920 
1921  // Loop over all edges of element i
1922  for(j = 0; j < exp->GetNedges(); ++j)
1923  {
1924  meshEdgeId = exp->GetGeom()->GetEid(j);
1925  pIt = perEdges.find(meshEdgeId);
1926  if (pIt != perEdges.end())
1927  {
1928  for (k = 0; k < pIt->second.size(); ++k)
1929  {
1930  meshEdgeId = min(meshEdgeId, pIt->second[k].id);
1931  }
1932  }
1933 
1934  edgeOrient = exp->GetGeom()->GetEorient(j);
1935  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1936  dof = exp->GetEdgeNcoeffs(j)-2;
1937 
1938  // Set the global DOF's for the interior modes of edge j
1939  for(k = 0; k < dof; ++k)
1940  {
1941  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
1942  m_globalToUniversalMap[vGlobalId]
1943  = nVert + meshEdgeId * maxEdgeDof + k + 1;
1944  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
1945  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
1946  }
1947  }
1948 
1949  // Loop over all faces of element i
1950  for(j = 0; j < exp->GetNfaces(); ++j)
1951  {
1952  faceOrient = exp->GetGeom()->GetForient(j);
1953 
1954  meshFaceId = exp->GetGeom()->GetFid(j);
1955 
1956  pIt = perFaces.find(meshFaceId);
1957  if (pIt != perFaces.end())
1958  {
1959  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
1960  {
1961  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1962  }
1963  meshFaceId = min(meshFaceId, pIt->second[0].id);
1964  }
1965 
1966 
1967  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1968  dof = exp->GetFaceIntNcoeffs(j);
1969 
1970 
1971  for(k = 0; k < dof; ++k)
1972  {
1973  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
1974  m_globalToUniversalMap[vGlobalId]
1975  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
1976  + k + 1;
1977  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
1978 
1979  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
1980  }
1981  }
1982 
1983  // Add interior DOFs to complete universal numbering
1984  exp->GetInteriorMap(interiorMap);
1985  dof = interiorMap.num_elements();
1986  elementId = (exp->GetGeom())->GetGlobalID();
1987  for (k = 0; k < dof; ++k)
1988  {
1989  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
1990  m_globalToUniversalMap[vGlobalId]
1991  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
1992  }
1993  }
1994 
1995  // Set up the GSLib universal assemble mapping
1996  // Internal DOF do not participate in any data
1997  // exchange, so we keep these set to the special GSLib id=0 so
1998  // they are ignored.
1999  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2000  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2001  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2002  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2003  {
2004  tmp[i] = m_globalToUniversalMap[i];
2005  }
2006 
2007  m_gsh = Gs::Init(tmp, vCommRow);
2008  m_bndGsh = Gs::Init(tmp2, vCommRow);
2009  Gs::Unique(tmp, vCommRow);
2010  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2011  {
2012  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2013  }
2014  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2015  {
2016  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2017  }
2018  }
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
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 2322 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().

2325  {
2327  if(global.data() == loc.data())
2328  {
2329  local = Array<OneD, NekDouble>(local.num_elements(),local.data());
2330  }
2331  else
2332  {
2333  local = loc; // create reference
2334  }
2335  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2336 
2337  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2338 
2339  if(m_signChange)
2340  {
2341  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2342  }
2343  else
2344  {
2345  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2346  }
2347  UniversalAssemble(global);
2348  }
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 2350 of file AssemblyMapCG.cpp.

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

2353  {
2354  Assemble(loc.GetPtr(),global.GetPtr());
2355  }
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 2419 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2420  {
2421  return m_extraDirEdges;
2422  }
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 2379 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2380  {
2381  return m_fullSystemBandWidth;
2382  }
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 2209 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2210  {
2211  return m_globalToUniversalMap[i];
2212  }
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 2226 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2227  {
2228  return m_globalToUniversalMap;
2229  }
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 2214 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2215  {
2216  return m_globalToUniversalMapUnique[i];
2217  }
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 2232 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2233  {
2235  }
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 2204 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2205  {
2206  return m_localToGlobalMap[i];
2207  }
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 2220 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2221  {
2222  return m_localToGlobalMap;
2223  }
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 2237 of file AssemblyMapCG.cpp.

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

2239  {
2240  if(m_signChange)
2241  {
2242  return m_localToGlobalSign[i];
2243  }
2244  else
2245  {
2246  return 1.0;
2247  }
2248  }
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 2250 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2251  {
2252  return m_localToGlobalSign;
2253  }
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 2399 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2400  {
2401  return m_numDirEdges;
2402  }
int m_numDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2404 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2405  {
2406  return m_numDirFaces;
2407  }
int m_numDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2389 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2390  {
2391  return m_numNonDirEdgeModes;
2392  }
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 2409 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2410  {
2411  return m_numNonDirEdges;
2412  }
int m_numNonDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2394 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2395  {
2396  return m_numNonDirFaceModes;
2397  }
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 2414 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2415  {
2416  return m_numNonDirFaces;
2417  }
int m_numNonDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2384 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2385  {
2386  return m_numNonDirVertexModes;
2387  }
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 2290 of file AssemblyMapCG.cpp.

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

2293  {
2295  if(global.data() == loc.data())
2296  {
2297  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2298  }
2299  else
2300  {
2301  glo = global; // create reference
2302  }
2303 
2304 
2305  if(m_signChange)
2306  {
2307  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2308  }
2309  else
2310  {
2311  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2312  }
2313  }
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 2315 of file AssemblyMapCG.cpp.

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

2318  {
2319  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2320  }
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 2028 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().

2030  {
2031  AssemblyMapCGSharedPtr returnval;
2032 
2033  int i, j;
2034  int nverts = 0;
2035  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2036  = locexp.GetExp();
2037  int nelmts = exp->size();
2038 
2039  // Get Default Map and turn off any searched values.
2040  returnval = MemoryManager<AssemblyMapCG>
2042  returnval->m_solnType = solnType;
2043  returnval->m_preconType = eNull;
2044  returnval->m_maxStaticCondLevel = 0;
2045  returnval->m_signChange = false;
2046  returnval->m_comm = m_comm;
2047 
2048  // Count the number of vertices
2049  for (i = 0; i < nelmts; ++i)
2050  {
2051  nverts += (*exp)[i]->GetNverts();
2052  }
2053 
2054  returnval->m_numLocalCoeffs = nverts;
2055  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2056 
2057  // Store original global ids in this map
2058  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2059 
2060  int cnt = 0;
2061  int cnt1 = 0;
2062  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2063 
2064  // Set up local to global map;
2065  for (i = 0; i < nelmts; ++i)
2066  {
2067  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2068  {
2069  returnval->m_localToGlobalMap[cnt] =
2070  returnval->m_localToGlobalBndMap[cnt] =
2071  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2072  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2073 
2074  // Set up numLocalDirBndCoeffs
2075  if ((returnval->m_localToGlobalMap[cnt]) <
2077  {
2078  returnval->m_numLocalDirBndCoeffs++;
2079  }
2080  cnt++;
2081  }
2082  cnt1 += (*exp)[i]->GetNcoeffs();
2083  }
2084 
2085  cnt = 0;
2086  // Reset global numbering and count number of dofs
2087  for (i = 0; i < m_numGlobalCoeffs; ++i)
2088  {
2089  if (GlobCoeffs[i] != -1)
2090  {
2091  GlobCoeffs[i] = cnt++;
2092  }
2093  }
2094 
2095  // Set up number of globalCoeffs;
2096  returnval->m_numGlobalCoeffs = cnt;
2097 
2098  // Set up number of global Dirichlet boundary coefficients
2099  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2100  {
2101  if (GlobCoeffs[i] != -1)
2102  {
2103  returnval->m_numGlobalDirBndCoeffs++;
2104  }
2105  }
2106 
2107  // Set up global to universal map
2108  if (m_globalToUniversalMap.num_elements())
2109  {
2111  = m_session->GetComm()->GetRowComm();
2112  int nglocoeffs = returnval->m_numGlobalCoeffs;
2113  returnval->m_globalToUniversalMap
2114  = Array<OneD, int> (nglocoeffs);
2115  returnval->m_globalToUniversalMapUnique
2116  = Array<OneD, int> (nglocoeffs);
2117 
2118  // Reset local to global map and setup universal map
2119  for (i = 0; i < nverts; ++i)
2120  {
2121  cnt = returnval->m_localToGlobalMap[i];
2122  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2123 
2124  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2126  }
2127 
2128  Nektar::Array<OneD, long> tmp(nglocoeffs);
2129  Vmath::Zero(nglocoeffs, tmp, 1);
2130  for (unsigned int i = 0; i < nglocoeffs; ++i)
2131  {
2132  tmp[i] = returnval->m_globalToUniversalMap[i];
2133  }
2134  returnval->m_gsh = Gs::Init(tmp, vCommRow);
2135  Gs::Unique(tmp, vCommRow);
2136  for (unsigned int i = 0; i < nglocoeffs; ++i)
2137  {
2138  returnval->m_globalToUniversalMapUnique[i]
2139  = (tmp[i] >= 0 ? 1 : 0);
2140  }
2141  }
2142  else // not sure this option is ever needed.
2143  {
2144  for (i = 0; i < nverts; ++i)
2145  {
2146  cnt = returnval->m_localToGlobalMap[i];
2147  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2148  }
2149  }
2150 
2151  return returnval;
2152  }
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 2255 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().

2258  {
2260  if(global.data() == loc.data())
2261  {
2262  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2263  }
2264  else
2265  {
2266  local = loc; // create reference
2267  }
2268 
2269 
2270  if(m_signChange)
2271  {
2272  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2273  }
2274  else
2275  {
2276  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2277  }
2278 
2279  // ensure all values are unique by calling a max
2280  Gs::Gather(global, Gs::gs_max, m_gsh);
2281  }
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 2283 of file AssemblyMapCG.cpp.

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

2286  {
2287  LocalToGlobal(loc.GetPtr(),global.GetPtr());
2288  }
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 2357 of file AssemblyMapCG.cpp.

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

2359  {
2360  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2361  }
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 2363 of file AssemblyMapCG.cpp.

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

2365  {
2366  UniversalAssemble(pGlobal.GetPtr());
2367  }
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 2369 of file AssemblyMapCG.cpp.

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

2372  {
2373  Array<OneD, NekDouble> tmp(offset);
2374  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2375  UniversalAssemble(pGlobal);
2376  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2377  }
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:1016

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 141 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 134 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 112 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 114 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 116 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 138 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG().

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

Number of Dirichlet edges.

Definition at line 124 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirEdges().

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

Number of Dirichlet faces.

Definition at line 126 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 132 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 136 of file AssemblyMapCG.h.

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

Number of non Dirichlet edge modes.

Definition at line 120 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().

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

Number of Dirichlet edges.

Definition at line 128 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 122 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().

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

Number of Dirichlet faces.

Definition at line 130 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 118 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().