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

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

#include <AssemblyMapCG.h>

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

Public Member Functions

 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Default constructor. More...
 
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const BndCondExp &bndCondExp=NullExpListSharedPtrArray, const BndCond &bndConditions=SpatialDomains::NullBoundaryConditionShPtrArray, const bool checkIfSingular=false, const std::string variable="defaultVar", const PeriodicMap &periodicVerts=NullPeriodicMap, const PeriodicMap &periodicEdges=NullPeriodicMap, const PeriodicMap &periodicFaces=NullPeriodicMap)
 General constructor for expansions of all dimensions without boundary conditions. More...
 
virtual ~AssemblyMapCG ()
 Destructor. More...
 
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 GetMaxIterations () const
 
int GetSuccessiveRHS () const
 
int GetLowestStaticCondLevel () const
 

Protected Member Functions

int CreateGraph (const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, 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...
 
int m_maxIterations
 Maximum iterations for iterative solver. More...
 
NekDouble m_iterativeTolerance
 Tolerance for iterative solver. More...
 
int m_successiveRHS
 sucessive RHS for iterative solver More...
 
Gs::gs_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
ExpListSharedPtr
BndCondExp
 
typedef Array< OneD, const
SpatialDomains::BoundaryConditionShPtr
BndCond
 

Detailed Description

Constructs mappings for the C0 scalar continuous Galerkin formulation.

Mappings are created for three possible global solution types:

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

Definition at line 71 of file AssemblyMapCG.h.

Member Typedef Documentation

Definition at line 75 of file AssemblyMapCG.h.

Definition at line 73 of file AssemblyMapCG.h.

Constructor & Destructor Documentation

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

Default constructor.

Definition at line 69 of file AssemblyMapCG.cpp.

References m_maxStaticCondLevel.

71  :
72  AssemblyMap(pSession,variable)
73  {
74  pSession->LoadParameter(
75  "MaxStaticCondLevel", m_maxStaticCondLevel, 100);
76  }
int m_maxStaticCondLevel
Maximum static condensation level.
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:77
Nektar::MultiRegions::AssemblyMapCG::AssemblyMapCG ( const LibUtilities::SessionReaderSharedPtr pSession,
const int  numLocalCoeffs,
const ExpList locExp,
const BndCondExp bndCondExp = NullExpListSharedPtrArray,
const BndCond bndConditions = SpatialDomains::NullBoundaryConditionShPtrArray,
const bool  checkIfSingular = false,
const std::string  variable = "defaultVar",
const PeriodicMap periodicVerts = NullPeriodicMap,
const PeriodicMap periodicEdges = NullPeriodicMap,
const PeriodicMap periodicFaces = NullPeriodicMap 
)

General constructor for expansions of all dimensions without boundary conditions.

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

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

Definition at line 1259 of file AssemblyMapCG.cpp.

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

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

Destructor.

Definition at line 1848 of file AssemblyMapCG.cpp.

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

1849  {
1852  }
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:205

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

2310  {
2311  int i,j;
2312  int cnt = 0;
2313  int locSize;
2314  int maxId;
2315  int minId;
2316  int bwidth = -1;
2317  for(i = 0; i < m_numPatches; ++i)
2318  {
2320  maxId = -1;
2321  minId = m_numLocalCoeffs+1;
2322  for(j = 0; j < locSize; j++)
2323  {
2325  {
2326  if(m_localToGlobalMap[cnt+j] > maxId)
2327  {
2328  maxId = m_localToGlobalMap[cnt+j];
2329  }
2330 
2331  if(m_localToGlobalMap[cnt+j] < minId)
2332  {
2333  minId = m_localToGlobalMap[cnt+j];
2334  }
2335  }
2336  }
2337  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2338 
2339  cnt+=locSize;
2340  }
2341 
2342  m_fullSystemBandWidth = bwidth;
2343  }
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
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:388
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:390
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:386
int Nektar::MultiRegions::AssemblyMapCG::CreateGraph ( const ExpList locExp,
const BndCondExp bndCondExp,
const Array< OneD, const BndCond > &  bndConditions,
const bool  checkIfSystemSingular,
const PeriodicMap periodicVerts,
const PeriodicMap periodicEdges,
const PeriodicMap periodicFaces,
DofGraph graph,
BottomUpSubStructuredGraphSharedPtr bottomUpGraph,
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::Finalise(), Gs::Gather(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), Nektar::MultiRegions::GlobalSysSolnTypeMap, Gs::gs_add, Vmath::Imax(), Gs::Init(), Nektar::iterator, Nektar::MultiRegions::AssemblyMap::m_comm, m_extraDirEdges, Nektar::MultiRegions::AssemblyMap::m_lowestStaticCondLevel, m_numDirEdges, m_numDirFaces, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs, m_numLocalBndCondCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs, m_numNonDirEdgeModes, m_numNonDirEdges, m_numNonDirFaceModes, m_numNonDirFaces, m_numNonDirVertexModes, Nektar::MultiRegions::AssemblyMap::m_session, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_systemSingular, Nektar::MultiRegions::MultiLevelBisectionReordering(), Nektar::MultiRegions::NoReordering(), Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceMin, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

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

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  if (EdgeSize.count(meshEdgeId) > 0)
588  {
589  EdgeSize[meshEdgeId] =
590  min(EdgeSize[meshEdgeId],
591  exp->GetEdgeNcoeffs(j) - 2);
592  }
593  else
594  {
595  EdgeSize[meshEdgeId] = exp->GetEdgeNcoeffs(j) - 2;
596  }
597  }
598 
599  nFaces = exp->GetNfaces();
600  faceCnt = 0;
601  for(j = 0; j < nFaces; ++j)
602  {
603  meshFaceId = exp->GetGeom()->GetFid(j);
604  if (FaceSize.count(meshFaceId) > 0)
605  {
606  FaceSize[meshFaceId] =
607  min(FaceSize[meshFaceId],
608  exp->GetFaceIntNcoeffs(j));
609  }
610  else
611  {
612  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
613  }
614  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
615  }
616  }
617 
618  /// - Periodic vertices
619  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
620  {
621  meshVertId = pIt->first;
622 
623  // This periodic vertex is joined to a Dirichlet condition.
624  if (graph[0].count(pIt->first) != 0)
625  {
626  for (i = 0; i < pIt->second.size(); ++i)
627  {
628  meshVertId2 = pIt->second[i].id;
629  if (graph[0].count(meshVertId2) == 0 &&
630  pIt->second[i].isLocal)
631  {
632  graph[0][meshVertId2] =
633  graph[0][meshVertId];
634  }
635  }
636  continue;
637  }
638 
639  // One of the attached vertices is Dirichlet.
640  bool isDirichlet = false;
641  for (i = 0; i < pIt->second.size(); ++i)
642  {
643  if (!pIt->second[i].isLocal)
644  {
645  continue;
646  }
647 
648  meshVertId2 = pIt->second[i].id;
649  if (graph[0].count(meshVertId2) > 0)
650  {
651  isDirichlet = true;
652  break;
653  }
654  }
655 
656  if (isDirichlet)
657  {
658  graph[0][meshVertId] =
659  graph[0][pIt->second[i].id];
660 
661  for (j = 0; j < pIt->second.size(); ++j)
662  {
663  meshVertId2 = pIt->second[i].id;
664  if (j == i || !pIt->second[j].isLocal ||
665  graph[0].count(meshVertId2) > 0)
666  {
667  continue;
668  }
669 
670  graph[0][meshVertId2] =
671  graph[0][pIt->second[i].id];
672  }
673 
674  continue;
675  }
676 
677  // Otherwise, see if a vertex ID has already been set.
678  for (i = 0; i < pIt->second.size(); ++i)
679  {
680  if (!pIt->second[i].isLocal)
681  {
682  continue;
683  }
684 
685  if (tempGraph[0].count(pIt->second[i].id) > 0)
686  {
687  break;
688  }
689  }
690 
691  if (i == pIt->second.size())
692  {
693  tempGraph[0][meshVertId] = tempGraphVertId++;
695  }
696  else
697  {
698  tempGraph[0][meshVertId] = tempGraph[0][pIt->second[i].id];
699  }
700  }
701 
702  // Store the temporary graph vertex id's of all element edges and
703  // vertices in these 3 arrays below
704  localVerts = Array<OneD, int>(nTotalVerts,-1);
705  localEdges = Array<OneD, int>(nTotalEdges,-1);
706  localFaces = Array<OneD, int>(nTotalFaces,-1);
707 
708  // Set up vertex numbering
709  for(i = 0; i < locExpVector.size(); ++i)
710  {
711  exp = locExpVector[i];
712  vertCnt = 0;
713  nVerts = exp->GetNverts();
714  for(j = 0; j < nVerts; ++j)
715  {
716  meshVertId = exp->GetGeom()->GetVid(j);
717  if(graph[0].count(meshVertId) == 0)
718  {
719  if(tempGraph[0].count(meshVertId) == 0)
720  {
721  boost::add_vertex(boostGraphObj);
722  tempGraph[0][meshVertId] = tempGraphVertId++;
724  }
725  localVerts[localVertOffset+vertCnt++] = tempGraph[0][meshVertId];
726  vwgts_map[ tempGraph[0][meshVertId] ] = 1;
727  }
728  }
729 
730  localVertOffset+=nVerts;
731  }
732 
733  /// - Periodic edges
734  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
735  {
736  meshEdgeId = pIt->first;
737 
738  // This periodic edge is joined to a Dirichlet condition.
739  if (graph[1].count(pIt->first) != 0)
740  {
741  for (i = 0; i < pIt->second.size(); ++i)
742  {
743  meshEdgeId2 = pIt->second[i].id;
744  if (graph[1].count(meshEdgeId2) == 0 &&
745  pIt->second[i].isLocal)
746  {
747  graph[1][meshEdgeId2] =
748  graph[1][meshEdgeId];
749  }
750  }
751  continue;
752  }
753 
754  // One of the attached edges is Dirichlet.
755  bool isDirichlet = false;
756  for (i = 0; i < pIt->second.size(); ++i)
757  {
758  if (!pIt->second[i].isLocal)
759  {
760  continue;
761  }
762 
763  meshEdgeId2 = pIt->second[i].id;
764  if (graph[1].count(meshEdgeId2) > 0)
765  {
766  isDirichlet = true;
767  break;
768  }
769  }
770 
771  if (isDirichlet)
772  {
773  graph[1][meshEdgeId] =
774  graph[1][pIt->second[i].id];
775 
776  for (j = 0; j < pIt->second.size(); ++j)
777  {
778  meshEdgeId2 = pIt->second[i].id;
779  if (j == i || !pIt->second[j].isLocal ||
780  graph[1].count(meshEdgeId2) > 0)
781  {
782  continue;
783  }
784 
785  graph[1][meshEdgeId2] =
786  graph[1][pIt->second[i].id];
787  }
788 
789  continue;
790  }
791 
792  // Otherwise, see if a edge ID has already been set.
793  for (i = 0; i < pIt->second.size(); ++i)
794  {
795  if (!pIt->second[i].isLocal)
796  {
797  continue;
798  }
799 
800  if (tempGraph[1].count(pIt->second[i].id) > 0)
801  {
802  break;
803  }
804  }
805 
806  if (i == pIt->second.size())
807  {
808  tempGraph[1][meshEdgeId] = tempGraphVertId++;
809  m_numNonDirEdgeModes += EdgeSize[meshEdgeId];
811  }
812  else
813  {
814  tempGraph[1][meshEdgeId] = tempGraph[1][pIt->second[i].id];
815  }
816  }
817 
818  int nEdgeIntCoeffs, nFaceIntCoeffs;
819 
820  // Set up edge numbering
821  for(i = 0; i < locExpVector.size(); ++i)
822  {
823  exp = locExpVector[i];
824  edgeCnt = 0;
825  nEdges = exp->GetNedges();
826 
827  for(j = 0; j < nEdges; ++j)
828  {
829  nEdgeIntCoeffs = exp->GetEdgeNcoeffs(j) - 2;
830  meshEdgeId = exp->GetGeom()->GetEid(j);
831  if(graph[1].count(meshEdgeId) == 0)
832  {
833  if(tempGraph[1].count(meshEdgeId) == 0)
834  {
835  boost::add_vertex(boostGraphObj);
836  tempGraph[1][meshEdgeId] = tempGraphVertId++;
837  m_numNonDirEdgeModes+=nEdgeIntCoeffs;
838 
840  }
841  localEdges[localEdgeOffset+edgeCnt++] = tempGraph[1][meshEdgeId];
842  vwgts_map[ tempGraph[1][meshEdgeId] ] = nEdgeIntCoeffs;
843  }
844  }
845 
846  localEdgeOffset+=nEdges;
847  }
848 
849  /// - Periodic faces
850  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
851  {
852  if (!pIt->second[0].isLocal)
853  {
854  // The face mapped to is on another process.
855  meshFaceId = pIt->first;
856  ASSERTL0(graph[2].count(meshFaceId) == 0,
857  "This periodic boundary edge has been specified before");
858  tempGraph[2][meshFaceId] = tempGraphVertId++;
859  nFaceIntCoeffs = FaceSize[meshFaceId];
860  m_numNonDirFaceModes+=nFaceIntCoeffs;
862  }
863  else if (pIt->first < pIt->second[0].id)
864  {
865  ASSERTL0(graph[2].count(pIt->first) == 0,
866  "This periodic boundary face has been specified before");
867  ASSERTL0(graph[2].count(pIt->second[0].id) == 0,
868  "This periodic boundary face has been specified before");
869 
870  tempGraph[2][pIt->first] = tempGraphVertId;
871  tempGraph[2][pIt->second[0].id] = tempGraphVertId++;
872  nFaceIntCoeffs = FaceSize[pIt->first];
873  m_numNonDirFaceModes+=nFaceIntCoeffs;
875  }
876  }
877 
878  // setup face numbering
879  for(i = 0; i < locExpVector.size(); ++i)
880  {
881  exp = locExpVector[i];
882  nFaces = exp->GetNfaces();
883  faceCnt = 0;
884  for(j = 0; j < nFaces; ++j)
885  {
886  nFaceIntCoeffs = exp->GetFaceIntNcoeffs(j);
887  meshFaceId = exp->GetGeom()->GetFid(j);
888  if(graph[2].count(meshFaceId) == 0)
889  {
890  if(tempGraph[2].count(meshFaceId) == 0)
891  {
892  boost::add_vertex(boostGraphObj);
893  tempGraph[2][meshFaceId] = tempGraphVertId++;
894  m_numNonDirFaceModes+=nFaceIntCoeffs;
895 
897  }
898  localFaces[localFaceOffset+faceCnt++] = tempGraph[2][meshFaceId];
899  vwgts_map[ tempGraph[2][meshFaceId] ] = nFaceIntCoeffs;
900  }
901  }
902  m_numLocalBndCoeffs += exp->NumBndryCoeffs();
903 
904  localFaceOffset+=nFaces;
905  }
906 
907  localVertOffset=0;
908  localEdgeOffset=0;
909  localFaceOffset=0;
910  for(i = 0; i < locExpVector.size(); ++i)
911  {
912  exp = locExpVector[i];
913  nVerts = exp->GetNverts();
914  nEdges = exp->GetNedges();
915  nFaces = exp->GetNfaces();
916 
917  // Now loop over all local faces, edges and vertices of this
918  // element and define that all other faces, edges and verices of
919  // this element are adjacent to them.
920 
921  // Vertices
922  for(j = 0; j < nVerts; j++)
923  {
924  if(localVerts[j+localVertOffset]==-1)
925  {
926  break;
927  }
928  // associate to other vertices
929  for(k = 0; k < nVerts; k++)
930  {
931  if(localVerts[k+localVertOffset]==-1)
932  {
933  break;
934  }
935  if(k!=j)
936  {
937  boost::add_edge( (size_t) localVerts[j+localVertOffset],
938  (size_t) localVerts[k+localVertOffset],boostGraphObj);
939  }
940  }
941  // associate to other edges
942  for(k = 0; k < nEdges; k++)
943  {
944  if(localEdges[k+localEdgeOffset]==-1)
945  {
946  break;
947  }
948  boost::add_edge( (size_t) localVerts[j+localVertOffset],
949  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
950  }
951  // associate to other faces
952  for(k = 0; k < nFaces; k++)
953  {
954  if(localFaces[k+localFaceOffset]==-1)
955  {
956  break;
957  }
958  boost::add_edge( (size_t) localVerts[j+localVertOffset],
959  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
960  }
961  }
962 
963  // Edges
964  for(j = 0; j < nEdges; j++)
965  {
966  if(localEdges[j+localEdgeOffset]==-1)
967  {
968  break;
969  }
970  // Associate to other edges
971  for(k = 0; k < nEdges; k++)
972  {
973  if(localEdges[k+localEdgeOffset]==-1)
974  {
975  break;
976  }
977  if(k!=j)
978  {
979  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
980  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
981  }
982  }
983  // Associate to vertices
984  for(k = 0; k < nVerts; k++)
985  {
986  if(localVerts[k+localVertOffset]==-1)
987  {
988  break;
989  }
990  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
991  (size_t) localVerts[k+localVertOffset],boostGraphObj);
992  }
993  // Associate to faces
994  for(k = 0; k < nFaces; k++)
995  {
996  if(localFaces[k+localFaceOffset]==-1)
997  {
998  break;
999  }
1000  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1001  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1002  }
1003  }
1004 
1005  // Faces
1006  for(j = 0; j < nFaces; j++)
1007  {
1008  if(localFaces[j+localFaceOffset]==-1)
1009  {
1010  break;
1011  }
1012  // Associate to other faces
1013  for(k = 0; k < nFaces; k++)
1014  {
1015  if(localFaces[k+localFaceOffset]==-1)
1016  {
1017  break;
1018  }
1019  if(k!=j)
1020  {
1021  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1022  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1023  }
1024  }
1025  // Associate to vertices
1026  for(k = 0; k < nVerts; k++)
1027  {
1028  if(localVerts[k+localVertOffset]==-1)
1029  {
1030  break;
1031  }
1032  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1033  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1034  }
1035  // Associate to edges
1036  for(k = 0; k < nEdges; k++)
1037  {
1038  if(localEdges[k+localEdgeOffset]==-1)
1039  {
1040  break;
1041  }
1042  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1043  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1044  }
1045  }
1046 
1047  localVertOffset+=nVerts;
1048  localEdgeOffset+=nEdges;
1049  localFaceOffset+=nFaces;
1050  }
1051 
1052  // Container to store vertices of the graph which correspond to
1053  // degrees of freedom along the boundary.
1054  set<int> partVerts;
1055 
1058  {
1059  vector<long> procVerts, procEdges, procFaces;
1060  set <int> foundVerts, foundEdges, foundFaces;
1061 
1062  // Loop over element and construct the procVerts and procEdges
1063  // vectors, which store the geometry IDs of mesh vertices and
1064  // edges respectively which are local to this process.
1065  for(i = cnt = 0; i < locExpVector.size(); ++i)
1066  {
1067  int elmtid = locExp.GetOffset_Elmt_Id(i);
1068  exp = locExpVector[elmtid];
1069  for (j = 0; j < exp->GetNverts(); ++j)
1070  {
1071  int vid = exp->GetGeom()->GetVid(j)+1;
1072  if (foundVerts.count(vid) == 0)
1073  {
1074  procVerts.push_back(vid);
1075  foundVerts.insert(vid);
1076  }
1077  }
1078 
1079  for (j = 0; j < exp->GetNedges(); ++j)
1080  {
1081  int eid = exp->GetGeom()->GetEid(j)+1;
1082 
1083  if (foundEdges.count(eid) == 0)
1084  {
1085  procEdges.push_back(eid);
1086  foundEdges.insert(eid);
1087  }
1088  }
1089 
1090  for (j = 0; j < exp->GetNfaces(); ++j)
1091  {
1092  int fid = exp->GetGeom()->GetFid(j)+1;
1093 
1094  if (foundFaces.count(fid) == 0)
1095  {
1096  procFaces.push_back(fid);
1097  foundFaces.insert(fid);
1098  }
1099  }
1100  }
1101 
1102  int unique_verts = foundVerts.size();
1103  int unique_edges = foundEdges.size();
1104  int unique_faces = foundFaces.size();
1105 
1106  // Now construct temporary GS objects. These will be used to
1107  // populate the arrays tmp3 and tmp4 with the multiplicity of
1108  // the vertices and edges respectively to identify those
1109  // vertices and edges which are located on partition boundary.
1110  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1111  Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm);
1112  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1113  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1114  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1115  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1116  Gs::Finalise(tmp1);
1117 
1118  if (unique_edges > 0)
1119  {
1120  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1121  Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
1122  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1123  Gs::Finalise(tmp2);
1124  }
1125 
1126  if (unique_faces > 0)
1127  {
1128  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1129  Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
1130  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1131  Gs::Finalise(tmp3);
1132  }
1133 
1134  // Finally, fill the partVerts set with all non-Dirichlet
1135  // vertices which lie on a partition boundary.
1136  for (i = 0; i < unique_verts; ++i)
1137  {
1138  if (tmp4[i] > 1.0)
1139  {
1140  if (graph[0].count(procVerts[i]-1) == 0)
1141  {
1142  partVerts.insert(tempGraph[0][procVerts[i]-1]);
1143  }
1144  }
1145  }
1146 
1147  for (i = 0; i < unique_edges; ++i)
1148  {
1149  if (tmp5[i] > 1.0)
1150  {
1151  if (graph[1].count(procEdges[i]-1) == 0)
1152  {
1153  partVerts.insert(tempGraph[1][procEdges[i]-1]);
1154  }
1155  }
1156  }
1157 
1158  for (i = 0; i < unique_faces; ++i)
1159  {
1160  if (tmp6[i] > 1.0)
1161  {
1162  if (graph[2].count(procFaces[i]-1) == 0)
1163  {
1164  partVerts.insert(tempGraph[2][procFaces[i]-1]);
1165  }
1166  }
1167  }
1168  }
1169 
1170  int nGraphVerts = tempGraphVertId;
1171  Array<OneD, int> perm(nGraphVerts);
1172  Array<OneD, int> iperm(nGraphVerts);
1173  Array<OneD, int> vwgts(nGraphVerts);
1174  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1175  for(i = 0; i < nGraphVerts; ++i)
1176  {
1177  vwgts[i] = vwgts_map[i];
1178  }
1179 
1180  if(nGraphVerts)
1181  {
1182  switch(m_solnType)
1183  {
1184  case eDirectFullMatrix:
1185  case eIterativeFull:
1186  case eIterativeStaticCond:
1187  case ePETScStaticCond:
1188  case ePETScFullMatrix:
1189  case eXxtFullMatrix:
1190  case eXxtStaticCond:
1191  {
1192  NoReordering(boostGraphObj,perm,iperm);
1193  break;
1194  }
1195 
1196  case eDirectStaticCond:
1197  {
1198  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1199  break;
1200  }
1201 
1206  {
1208  boostGraphObj, perm, iperm, bottomUpGraph,
1209  partVerts, mdswitch);
1210  break;
1211  }
1212  default:
1213  {
1214  ASSERTL0(false,
1215  "Unrecognised solution type: " + std::string(
1217  }
1218  }
1219  }
1220 
1221  // For parallel multi-level static condensation determine the lowest
1222  // static condensation level amongst processors.
1225  {
1226  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1227  vComm->AllReduce(m_lowestStaticCondLevel,
1229  }
1230  else
1231  {
1233  }
1234 
1235  /**
1236  * STEP 4: Fill the #graph[0] and
1237  * #graph[1] with the optimal ordering from boost.
1238  */
1240  for(mapIt = tempGraph[0].begin(); mapIt != tempGraph[0].end(); mapIt++)
1241  {
1242  graph[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1243  }
1244  for(mapIt = tempGraph[1].begin(); mapIt != tempGraph[1].end(); mapIt++)
1245  {
1246  graph[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1247  }
1248  for(mapIt = tempGraph[2].begin(); mapIt != tempGraph[2].end(); mapIt++)
1249  {
1250  graph[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1251  }
1252 
1253  return nGraphVerts;
1254  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:320
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:306
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:223
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:205
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:741
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:363
int m_numDirFaces
Number of Dirichlet faces.
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:397
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:316
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:312
int m_numDirEdges
Number of Dirichlet edges.
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:148
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:303
int m_numNonDirFaces
Number of Dirichlet faces.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
map<int, vector<ExtraDirDof> >& Nektar::MultiRegions::AssemblyMapCG::GetExtraDirDofs ( )
inline

Definition at line 108 of file AssemblyMapCG.h.

References m_extraDirDofs.

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

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

Definition at line 1964 of file AssemblyMapCG.cpp.

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

Referenced by AssemblyMapCG().

1969  {
1971  int nVert = 0;
1972  int nEdge = 0;
1973  int nFace = 0;
1974  int maxEdgeDof = 0;
1975  int maxFaceDof = 0;
1976  int maxIntDof = 0;
1977  int dof = 0;
1978  int cnt;
1979  int i,j,k;
1980  int meshVertId;
1981  int meshEdgeId;
1982  int meshFaceId;
1983  int elementId;
1984  int vGlobalId;
1985  int maxBndGlobalId = 0;
1986  StdRegions::Orientation edgeOrient;
1987  StdRegions::Orientation faceOrient;
1988  Array<OneD, unsigned int> edgeInteriorMap;
1989  Array<OneD, int> edgeInteriorSign;
1990  Array<OneD, unsigned int> faceInteriorMap;
1991  Array<OneD, int> faceInteriorSign;
1992  Array<OneD, unsigned int> interiorMap;
1993  PeriodicMap::const_iterator pIt;
1994 
1995  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1996  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
1997 
2002 
2003  // Loop over all the elements in the domain to gather mesh data
2004  for(i = 0; i < locExpVector.size(); ++i)
2005  {
2006  exp = locExpVector[i];
2007  nVert += exp->GetNverts();
2008  nEdge += exp->GetNedges();
2009  nFace += exp->GetNfaces();
2010  // Loop over all edges (and vertices) of element i
2011  for(j = 0; j < exp->GetNedges(); ++j)
2012  {
2013  dof = exp->GetEdgeNcoeffs(j)-2;
2014  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2015  }
2016  for(j = 0; j < exp->GetNfaces(); ++j)
2017  {
2018  dof = exp->GetFaceIntNcoeffs(j);
2019  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2020  }
2021  exp->GetInteriorMap(interiorMap);
2022  dof = interiorMap.num_elements();
2023  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2024  }
2025 
2026  // Tell other processes about how many dof we have
2027  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2028  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2029  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2030  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2031  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2032  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2033 
2034  // Assemble global to universal mapping for this process
2035  for(i = 0; i < locExpVector.size(); ++i)
2036  {
2037  exp = locExpVector[i];
2038  cnt = locExp.GetCoeff_Offset(i);
2039 
2040  // Loop over all vertices of element i
2041  for(j = 0; j < exp->GetNverts(); ++j)
2042  {
2043  meshVertId = exp->GetGeom()->GetVid(j);
2044  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2045 
2046  pIt = perVerts.find(meshVertId);
2047  if (pIt != perVerts.end())
2048  {
2049  for (k = 0; k < pIt->second.size(); ++k)
2050  {
2051  meshVertId = min(meshVertId, pIt->second[k].id);
2052  }
2053  }
2054 
2055  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2056  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2057  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2058  }
2059 
2060  // Loop over all edges of element i
2061  for(j = 0; j < exp->GetNedges(); ++j)
2062  {
2063  meshEdgeId = exp->GetGeom()->GetEid(j);
2064  pIt = perEdges.find(meshEdgeId);
2065  edgeOrient = exp->GetGeom()->GetEorient(j);
2066 
2067  if (pIt != perEdges.end())
2068  {
2069  pair<int, StdRegions::Orientation> idOrient =
2071  meshEdgeId, edgeOrient, pIt->second);
2072  meshEdgeId = idOrient.first;
2073  edgeOrient = idOrient.second;
2074  }
2075 
2076  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
2077  dof = exp->GetEdgeNcoeffs(j)-2;
2078 
2079  // Set the global DOF's for the interior modes of edge j
2080  // run backwards because of variable P case "ghost" modes
2081  for(k = dof-1; k >= 0; --k)
2082  {
2083  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2084  m_globalToUniversalMap[vGlobalId]
2085  = nVert + meshEdgeId * maxEdgeDof + k + 1;
2086  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2087  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2088  }
2089  }
2090 
2091  // Loop over all faces of element i
2092  for(j = 0; j < exp->GetNfaces(); ++j)
2093  {
2094  faceOrient = exp->GetGeom()->GetForient(j);
2095 
2096  meshFaceId = exp->GetGeom()->GetFid(j);
2097 
2098  pIt = perFaces.find(meshFaceId);
2099  if (pIt != perFaces.end())
2100  {
2101  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2102  {
2103  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
2104  }
2105  meshFaceId = min(meshFaceId, pIt->second[0].id);
2106  }
2107 
2108 
2109  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2110  dof = exp->GetFaceIntNcoeffs(j);
2111 
2112 
2113  for(k = dof-1; k >= 0; --k)
2114  {
2115  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2116  m_globalToUniversalMap[vGlobalId]
2117  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2118  + k + 1;
2119  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2120 
2121  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2122  }
2123  }
2124 
2125  // Add interior DOFs to complete universal numbering
2126  exp->GetInteriorMap(interiorMap);
2127  dof = interiorMap.num_elements();
2128  elementId = (exp->GetGeom())->GetGlobalID();
2129  for (k = 0; k < dof; ++k)
2130  {
2131  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2132  m_globalToUniversalMap[vGlobalId]
2133  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2134  }
2135  }
2136 
2137  // Set up the GSLib universal assemble mapping
2138  // Internal DOF do not participate in any data
2139  // exchange, so we keep these set to the special GSLib id=0 so
2140  // they are ignored.
2141  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2142  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2143  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2144  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2145  {
2146  tmp[i] = m_globalToUniversalMap[i];
2147  }
2148 
2149  m_gsh = Gs::Init(tmp, vCommRow);
2150  m_bndGsh = Gs::Init(tmp2, vCommRow);
2151  Gs::Unique(tmp, vCommRow);
2152  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2153  {
2154  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2155  }
2156  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2157  {
2158  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2159  }
2160  }
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
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:183
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:358
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:360
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
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 2464 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().

2467  {
2469  if(global.data() == loc.data())
2470  {
2471  local = Array<OneD, NekDouble>(local.num_elements(),local.data());
2472  }
2473  else
2474  {
2475  local = loc; // create reference
2476  }
2477  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2478 
2479  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2480 
2481  if(m_signChange)
2482  {
2483  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2484  }
2485  else
2486  {
2487  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2488  }
2489  UniversalAssemble(global);
2490  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
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:695
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
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 2492 of file AssemblyMapCG.cpp.

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

2495  {
2496  Assemble(loc.GetPtr(),global.GetPtr());
2497  }
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 2561 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2562  {
2563  return m_extraDirEdges;
2564  }
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 2521 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2522  {
2523  return m_fullSystemBandWidth;
2524  }
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 2351 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2352  {
2353  return m_globalToUniversalMap[i];
2354  }
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 2368 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2369  {
2370  return m_globalToUniversalMap;
2371  }
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 2356 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2357  {
2358  return m_globalToUniversalMapUnique[i];
2359  }
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 2374 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2375  {
2377  }
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 2346 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2347  {
2348  return m_localToGlobalMap[i];
2349  }
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 2362 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2363  {
2364  return m_localToGlobalMap;
2365  }
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 2379 of file AssemblyMapCG.cpp.

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

2381  {
2382  if(m_signChange)
2383  {
2384  return m_localToGlobalSign[i];
2385  }
2386  else
2387  {
2388  return 1.0;
2389  }
2390  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
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 2392 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2393  {
2394  return m_localToGlobalSign;
2395  }
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 2541 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2542  {
2543  return m_numDirEdges;
2544  }
int m_numDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2546 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2547  {
2548  return m_numDirFaces;
2549  }
int m_numDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2531 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2532  {
2533  return m_numNonDirEdgeModes;
2534  }
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 2551 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2552  {
2553  return m_numNonDirEdges;
2554  }
int m_numNonDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2536 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2537  {
2538  return m_numNonDirFaceModes;
2539  }
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 2556 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2557  {
2558  return m_numNonDirFaces;
2559  }
int m_numNonDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2526 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2527  {
2528  return m_numNonDirVertexModes;
2529  }
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 2432 of file AssemblyMapCG.cpp.

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

2435  {
2437  if(global.data() == loc.data())
2438  {
2439  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2440  }
2441  else
2442  {
2443  glo = global; // create reference
2444  }
2445 
2446 
2447  if(m_signChange)
2448  {
2449  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2450  }
2451  else
2452  {
2453  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2454  }
2455  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
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:331
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 2457 of file AssemblyMapCG.cpp.

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

2460  {
2461  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2462  }
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 2170 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().

2172  {
2173  AssemblyMapCGSharedPtr returnval;
2174 
2175  int i, j;
2176  int nverts = 0;
2177  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2178  = locexp.GetExp();
2179  int nelmts = exp->size();
2180 
2181  // Get Default Map and turn off any searched values.
2182  returnval = MemoryManager<AssemblyMapCG>
2184  returnval->m_solnType = solnType;
2185  returnval->m_preconType = eNull;
2186  returnval->m_maxStaticCondLevel = 0;
2187  returnval->m_signChange = false;
2188  returnval->m_comm = m_comm;
2189 
2190  // Count the number of vertices
2191  for (i = 0; i < nelmts; ++i)
2192  {
2193  nverts += (*exp)[i]->GetNverts();
2194  }
2195 
2196  returnval->m_numLocalCoeffs = nverts;
2197  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2198 
2199  // Store original global ids in this map
2200  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2201 
2202  int cnt = 0;
2203  int cnt1 = 0;
2204  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2205 
2206  // Set up local to global map;
2207  for (i = 0; i < nelmts; ++i)
2208  {
2209  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2210  {
2211  returnval->m_localToGlobalMap[cnt] =
2212  returnval->m_localToGlobalBndMap[cnt] =
2213  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2214  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2215 
2216  // Set up numLocalDirBndCoeffs
2217  if ((returnval->m_localToGlobalMap[cnt]) <
2219  {
2220  returnval->m_numLocalDirBndCoeffs++;
2221  }
2222  cnt++;
2223  }
2224  cnt1 += (*exp)[i]->GetNcoeffs();
2225  }
2226 
2227  cnt = 0;
2228  // Reset global numbering and count number of dofs
2229  for (i = 0; i < m_numGlobalCoeffs; ++i)
2230  {
2231  if (GlobCoeffs[i] != -1)
2232  {
2233  GlobCoeffs[i] = cnt++;
2234  }
2235  }
2236 
2237  // Set up number of globalCoeffs;
2238  returnval->m_numGlobalCoeffs = cnt;
2239 
2240  // Set up number of global Dirichlet boundary coefficients
2241  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2242  {
2243  if (GlobCoeffs[i] != -1)
2244  {
2245  returnval->m_numGlobalDirBndCoeffs++;
2246  }
2247  }
2248 
2249  // Set up global to universal map
2250  if (m_globalToUniversalMap.num_elements())
2251  {
2253  = m_session->GetComm()->GetRowComm();
2254  int nglocoeffs = returnval->m_numGlobalCoeffs;
2255  returnval->m_globalToUniversalMap
2256  = Array<OneD, int> (nglocoeffs);
2257  returnval->m_globalToUniversalMapUnique
2258  = Array<OneD, int> (nglocoeffs);
2259 
2260  // Reset local to global map and setup universal map
2261  for (i = 0; i < nverts; ++i)
2262  {
2263  cnt = returnval->m_localToGlobalMap[i];
2264  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2265 
2266  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2268  }
2269 
2270  Nektar::Array<OneD, long> tmp(nglocoeffs);
2271  Vmath::Zero(nglocoeffs, tmp, 1);
2272  for (unsigned int i = 0; i < nglocoeffs; ++i)
2273  {
2274  tmp[i] = returnval->m_globalToUniversalMap[i];
2275  }
2276  returnval->m_gsh = Gs::Init(tmp, vCommRow);
2277  Gs::Unique(tmp, vCommRow);
2278  for (unsigned int i = 0; i < nglocoeffs; ++i)
2279  {
2280  returnval->m_globalToUniversalMapUnique[i]
2281  = (tmp[i] >= 0 ? 1 : 0);
2282  }
2283  }
2284  else // not sure this option is ever needed.
2285  {
2286  for (i = 0; i < nverts; ++i)
2287  {
2288  cnt = returnval->m_localToGlobalMap[i];
2289  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2290  }
2291  }
2292 
2293  return returnval;
2294  }
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
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:318
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:183
No Solution type specified.
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:303
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
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 2397 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().

2400  {
2402  if(global.data() == loc.data())
2403  {
2404  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2405  }
2406  else
2407  {
2408  local = loc; // create reference
2409  }
2410 
2411 
2412  if(m_signChange)
2413  {
2414  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2415  }
2416  else
2417  {
2418  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2419  }
2420 
2421  // ensure all values are unique by calling a max
2422  Gs::Gather(global, Gs::gs_max, m_gsh);
2423  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
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:223
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
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 2425 of file AssemblyMapCG.cpp.

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

2428  {
2429  LocalToGlobal(loc.GetPtr(),global.GetPtr());
2430  }
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 2499 of file AssemblyMapCG.cpp.

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

2501  {
2502  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2503  }
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:223
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( NekVector< NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2505 of file AssemblyMapCG.cpp.

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

2507  {
2508  UniversalAssemble(pGlobal.GetPtr());
2509  }
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 2511 of file AssemblyMapCG.cpp.

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

2514  {
2515  Array<OneD, NekDouble> tmp(offset);
2516  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2517  UniversalAssemble(pGlobal);
2518  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2519  }
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:1047

Member Data Documentation

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

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

Definition at line 148 of file AssemblyMapCG.h.

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

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

Extra dirichlet edges in parallel.

Definition at line 141 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetExtraDirEdges().

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

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

Definition at line 119 of file AssemblyMapCG.h.

Referenced by CalculateFullSystemBandWidth(), and v_GetFullSystemBandWidth().

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

Integer map of process coeffs to universal space.

Definition at line 121 of file AssemblyMapCG.h.

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

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

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

Definition at line 123 of file AssemblyMapCG.h.

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

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

Maximum static condensation level.

Definition at line 145 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG().

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

Number of Dirichlet edges.

Definition at line 131 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirEdges().

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

Number of Dirichlet faces.

Definition at line 133 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirFaces().

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

Number of local boundary condition coefficients.

Definition at line 139 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and CreateGraph().

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

Number of local boundary condition degrees of freedom.

Definition at line 143 of file AssemblyMapCG.h.

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

Number of non Dirichlet edge modes.

Definition at line 127 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().

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

Number of Dirichlet edges.

Definition at line 135 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdges().

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

Number of non Dirichlet face modes.

Definition at line 129 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().

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

Number of Dirichlet faces.

Definition at line 137 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaces().

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

Number of non Dirichlet vertex modes.

Definition at line 125 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().