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

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

#include <AssemblyMapCG.h>

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

Public Member Functions

 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Default constructor. More...
 
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const BndCondExp &bndCondExp=NullExpListSharedPtrArray, const BndCond &bndConditions=SpatialDomains::NullBoundaryConditionShPtrArray, const bool checkIfSingular=false, const std::string variable="defaultVar", const PeriodicMap &periodicVerts=NullPeriodicMap, const PeriodicMap &periodicEdges=NullPeriodicMap, const PeriodicMap &periodicFaces=NullPeriodicMap)
 General constructor for expansions of all dimensions without boundary conditions. More...
 
virtual ~AssemblyMapCG ()
 Destructor. More...
 
std::map< int, std::vector
< ExtraDirDof > > & 
GetExtraDirDofs ()
 
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
 AssemblyMap ()
 Default constructor. More...
 
 AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Constructor with a communicator. More...
 
 AssemblyMap (AssemblyMap *oldLevelMap, const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
 Constructor for next level in multi-level static condensation. More...
 
virtual ~AssemblyMap ()
 Destructor. More...
 
LibUtilities::CommSharedPtr GetComm ()
 Retrieves the communicator. More...
 
size_t GetHash () const
 Retrieves the hash of this map. More...
 
int GetLocalToGlobalMap (const int i) const
 
int GetGlobalToUniversalMap (const int i) const
 
int GetGlobalToUniversalMapUnique (const int i) const
 
const Array< OneD, const int > & GetLocalToGlobalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMap ()
 
const Array< OneD, const int > & GetGlobalToUniversalMapUnique ()
 
NekDouble GetLocalToGlobalSign (const int i) const
 
const Array< OneD, NekDouble > & GetLocalToGlobalSign () const
 
void LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
 
void LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm=true) 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, bool printHeader=true) 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, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
 
void SetUpUniversalC0ContMap (const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
 
void CalculateFullSystemBandWidth ()
 Calculate the bandwith of the full matrix system. More...
 
virtual int v_GetLocalToGlobalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMap (const int i) const
 
virtual int v_GetGlobalToUniversalMapUnique (const int i) const
 
virtual const Array< OneD,
const int > & 
v_GetLocalToGlobalMap ()
 
virtual const Array< OneD,
const int > & 
v_GetGlobalToUniversalMap ()
 
virtual const Array< OneD,
const int > & 
v_GetGlobalToUniversalMapUnique ()
 
virtual NekDouble v_GetLocalToGlobalSign (const int i) const
 
virtual const Array< OneD,
NekDouble > & 
v_GetLocalToGlobalSign () const
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
 
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm) 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...
 
std::map< int, std::vector
< ExtraDirDof > > 
m_extraDirDofs
 Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. More...
 
- Protected Attributes inherited from Nektar::MultiRegions::AssemblyMap
LibUtilities::SessionReaderSharedPtr m_session
 Session object. More...
 
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
size_t m_hash
 Hash for map. More...
 
int m_numLocalBndCoeffs
 Number of local boundary coefficients. More...
 
int m_numGlobalBndCoeffs
 Total number of global boundary coefficients. More...
 
int m_numLocalDirBndCoeffs
 Number of Local Dirichlet Boundary Coefficients. More...
 
int m_numGlobalDirBndCoeffs
 Number of Global Dirichlet Boundary Coefficients. More...
 
bool m_systemSingular
 Flag indicating if the system is singular or not. More...
 
int m_numLocalCoeffs
 Total number of local coefficients. More...
 
int m_numGlobalCoeffs
 Total number of global coefficients. More...
 
bool m_signChange
 Flag indicating if modes require sign reversal. More...
 
Array< OneD, int > m_localToGlobalBndMap
 Integer map of local boundary coeffs to global space. More...
 
Array< OneD, 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 72 of file AssemblyMapCG.cpp.

References m_maxStaticCondLevel.

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

General constructor for expansions of all dimensions without boundary conditions.

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

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

Definition at line 1295 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::LibUtilities::eModified_C, Nektar::MultiRegions::ePETScMultiLevelStaticCond, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Gs::Finalise(), Gs::Gather(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::LibUtilities::StdTriData::getNumberOfBndCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfBndCoefficients(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), 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_session, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_staticCondLevel, class_topology::P, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceSum, SetUpUniversalC0ContMap(), Nektar::MultiRegions::AssemblyMap::UniversalAssembleBnd(), and Vmath::Vmax().

1306  : AssemblyMap(pSession, variable)
1307  {
1308  int i, j, k, l;
1309  int p, q, numModes0, numModes1;
1310  int cnt = 0;
1311  int intDofCnt;
1312  int meshVertId, meshEdgeId, meshEdgeId2, meshFaceId, meshFaceId2;
1313  int globalId;
1314  int nEdgeInteriorCoeffs;
1315  int firstNonDirGraphVertId;
1316  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1319  StdRegions::Orientation edgeOrient;
1320  StdRegions::Orientation faceOrient;
1321  Array<OneD, unsigned int> edgeInteriorMap;
1322  Array<OneD, int> edgeInteriorSign;
1323  Array<OneD, unsigned int> faceInteriorMap;
1324  Array<OneD, int> faceInteriorSign;
1325  PeriodicMap::const_iterator pIt;
1326 
1327  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1328 
1329  bool verbose = m_session->DefinesCmdLineArgument("verbose");
1330 
1331  m_signChange = false;
1332 
1333  // Stores vertex, edge and face reordered vertices.
1334  DofGraph graph(3);
1335  DofGraph dofs(3);
1336  vector<map<int, int> > faceModes(2);
1337  map<int, LibUtilities::ShapeType> faceType;
1338 
1339  set<int> extraDirVerts, extraDirEdges;
1341 
1342  // Construct list of number of degrees of freedom for each vertex,
1343  // edge and face.
1344  for (i = 0; i < locExpVector.size(); ++i)
1345  {
1346  exp = locExpVector[i];
1347 
1348  for(j = 0; j < exp->GetNverts(); ++j)
1349  {
1350  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1351  }
1352 
1353  for(j = 0; j < exp->GetNedges(); ++j)
1354  {
1355  if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1356  {
1357  if (dofs[1][exp->GetGeom()->GetEid(j)] !=
1358  exp->GetEdgeNcoeffs(j)-2)
1359  {
1360  ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
1361  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B) ||
1362  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_C),
1363  "CG with variable order only available with modal expansion");
1364  }
1365  dofs[1][exp->GetGeom()->GetEid(j)] =
1366  min(dofs[1][exp->GetGeom()->GetEid(j)],
1367  exp->GetEdgeNcoeffs(j)-2);
1368  }
1369  else
1370  {
1371  dofs[1][exp->GetGeom()->GetEid(j)] =
1372  exp->GetEdgeNcoeffs(j) - 2;
1373  }
1374  }
1375 
1376  for(j = 0; j < exp->GetNfaces(); ++j)
1377  {
1378  faceOrient = exp->GetGeom()->GetForient(j);
1379  meshFaceId = exp->GetGeom()->GetFid(j);
1380  exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
1381 
1382  if (faceModes[0].count(meshFaceId) > 0)
1383  {
1384  faceModes[0][meshFaceId] =
1385  min(faceModes[0][meshFaceId], numModes0);
1386 
1387  faceModes[1][meshFaceId] =
1388  min(faceModes[1][meshFaceId], numModes1);
1389  }
1390  else
1391  {
1392  faceModes[0][meshFaceId] = numModes0;
1393  faceModes[1][meshFaceId] = numModes1;
1394 
1395  // Get shape of this face
1397  geom = boost::dynamic_pointer_cast<SpatialDomains::
1398  Geometry3D> (exp->GetGeom());
1399  faceType[meshFaceId] =
1400  geom->GetFace(j)->GetShapeType();
1401  }
1402  }
1403  }
1404 
1405  // Add non-local periodic dofs to the map
1406  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1407  {
1408  for (i = 0; i < pIt->second.size(); ++i)
1409  {
1410  meshEdgeId2 = pIt->second[i].id;
1411  if (dofs[1].count(meshEdgeId2) == 0)
1412  {
1413  dofs[1][meshEdgeId2] = 1e6;
1414  }
1415  }
1416  }
1417  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1418  {
1419  for (i = 0; i < pIt->second.size(); ++i)
1420  {
1421  meshFaceId2 = pIt->second[i].id;
1422  if (faceModes[0].count(meshFaceId2) == 0)
1423  {
1424  faceModes[0][meshFaceId2] = 1e6;
1425  faceModes[1][meshFaceId2] = 1e6;
1426  }
1427  }
1428  }
1429 
1430  // Now use information from all partitions to determine
1431  // the correct size
1432  map<int, int>::iterator dofIt, dofIt2;
1433  // edges
1434  Array<OneD, long> edgeId (dofs[1].size());
1435  Array<OneD, NekDouble> edgeDof (dofs[1].size());
1436  for(dofIt = dofs[1].begin(), i=0; dofIt != dofs[1].end(); dofIt++, i++)
1437  {
1438  edgeId[i] = dofIt->first + 1;
1439  edgeDof[i] = (NekDouble) dofIt->second;
1440  }
1441  Gs::gs_data *tmp = Gs::Init(edgeId, vComm, verbose);
1442  Gs::Gather(edgeDof, Gs::gs_min, tmp);
1443  Gs::Finalise(tmp);
1444  for (i=0; i < dofs[1].size(); i++)
1445  {
1446  dofs[1][edgeId[i]-1] = (int) (edgeDof[i]+0.5);
1447  }
1448  // Periodic edges
1449  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1450  {
1451  meshEdgeId = pIt->first;
1452  for (i = 0; i < pIt->second.size(); ++i)
1453  {
1454  meshEdgeId2 = pIt->second[i].id;
1455  if (dofs[1][meshEdgeId2] < dofs[1][meshEdgeId])
1456  {
1457  dofs[1][meshEdgeId] = dofs[1][meshEdgeId2];
1458  }
1459  }
1460  }
1461  // faces
1462  Array<OneD, long> faceId (faceModes[0].size());
1463  Array<OneD, NekDouble> faceP (faceModes[0].size());
1464  Array<OneD, NekDouble> faceQ (faceModes[0].size());
1465  for(dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin(),i=0;
1466  dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
1467  {
1468  faceId[i] = dofIt->first+1;
1469  faceP[i] = (NekDouble) dofIt->second;
1470  faceQ[i] = (NekDouble) dofIt2->second;
1471  }
1472  Gs::gs_data *tmp2 = Gs::Init(faceId, vComm, verbose);
1473  Gs::Gather(faceP, Gs::gs_min, tmp2);
1474  Gs::Gather(faceQ, Gs::gs_min, tmp2);
1475  Gs::Finalise(tmp2);
1476  for (i=0; i < faceModes[0].size(); i++)
1477  {
1478  faceModes[0][faceId[i]-1] = (int) (faceP[i]+0.5);
1479  faceModes[1][faceId[i]-1] = (int) (faceQ[i]+0.5);
1480  }
1481  // Periodic faces
1482  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1483  {
1484  meshFaceId = pIt->first;
1485  for (i = 0; i < pIt->second.size(); ++i)
1486  {
1487  meshFaceId2 = pIt->second[i].id;
1488  if (faceModes[0][meshFaceId2] < faceModes[0][meshFaceId])
1489  {
1490  faceModes[0][meshFaceId] = faceModes[0][meshFaceId2];
1491  }
1492  if (faceModes[1][meshFaceId2] < faceModes[1][meshFaceId])
1493  {
1494  faceModes[1][meshFaceId] = faceModes[1][meshFaceId2];
1495  }
1496  }
1497  }
1498  // Calculate number of dof in each face
1499  int P, Q;
1500  for (i=0; i < faceModes[0].size(); i++)
1501  {
1502  P = faceModes[0][faceId[i]-1];
1503  Q = faceModes[1][faceId[i]-1];
1504  if (faceType[faceId[i]-1] == LibUtilities::eQuadrilateral)
1505  {
1506  // Quad face
1507  dofs[2][faceId[i]-1] =
1510  }
1511  else
1512  {
1513  // Tri face
1514  dofs[2][faceId[i]-1] =
1517  }
1518  }
1519 
1520  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1521 
1522  // Note that nExtraDirichlet is not used in the logic below; it just
1523  // needs to be set so that the coupled solver in
1524  // IncNavierStokesSolver can work.
1525  int nExtraDirichlet;
1526  int mdswitch;
1527  m_session->LoadParameter(
1528  "MDSwitch", mdswitch, 10);
1529  int nGraphVerts =
1530  CreateGraph(locExp, bndCondExp, bndCondVec,
1531  checkIfSystemSingular, periodicVerts, periodicEdges,
1532  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1533  extraDirEdges, firstNonDirGraphVertId,
1534  nExtraDirichlet, mdswitch);
1535 
1536  /*
1537  * Set up an array which contains the offset information of the
1538  * different graph vertices.
1539  *
1540  * This basically means to identify to how many global degrees of
1541  * freedom the individual graph vertices correspond. Obviously,
1542  * the graph vertices corresponding to the mesh-vertices account
1543  * for a single global DOF. However, the graph vertices
1544  * corresponding to the element edges correspond to N-2 global DOF
1545  * where N is equal to the number of boundary modes on this edge.
1546  */
1547  Array<OneD, int> graphVertOffset(
1548  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1549 
1550  graphVertOffset[0] = 0;
1551 
1552  for(i = 0; i < locExpVector.size(); ++i)
1553  {
1554  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1555 
1556  for(j = 0; j < exp->GetNverts(); ++j)
1557  {
1558  meshVertId = exp->GetGeom()->GetVid(j);
1559  graphVertOffset[graph[0][meshVertId]+1] = 1;
1560  }
1561 
1562  for(j = 0; j < exp->GetNedges(); ++j)
1563  {
1564  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1565  meshEdgeId = exp->GetGeom()->GetEid(j);
1566  graphVertOffset[graph[1][meshEdgeId]+1]
1567  = dofs[1][meshEdgeId];
1568 
1569  bType = exp->GetEdgeBasisType(j);
1570 
1571  // need a sign vector for modal expansions if nEdgeCoeffs
1572  // >=3 (not 4 because of variable order case)
1573  if(nEdgeInteriorCoeffs &&
1574  (bType == LibUtilities::eModified_A ||
1575  bType == LibUtilities::eModified_B))
1576  {
1577  m_signChange = true;
1578  }
1579  }
1580 
1581  for(j = 0; j < exp->GetNfaces(); ++j)
1582  {
1583  meshFaceId = exp->GetGeom()->GetFid(j);
1584  graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
1585  }
1586  }
1587 
1588  for(i = 1; i < graphVertOffset.num_elements(); i++)
1589  {
1590  graphVertOffset[i] += graphVertOffset[i-1];
1591  }
1592 
1593  // Allocate the proper amount of space for the class-data
1594  m_numLocalCoeffs = numLocalCoeffs;
1595  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1596  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
1597  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
1599  // If required, set up the sign-vector
1600  if(m_signChange)
1601  {
1602  m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
1603  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
1604  m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD,NekDouble>(m_numLocalBndCondCoeffs,1.0);
1605  }
1606 
1607  m_staticCondLevel = 0;
1608  m_numPatches = locExpVector.size();
1609  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1610  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1611  for(i = 0; i < m_numPatches; ++i)
1612  {
1613  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1614  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1615  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1616  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1617  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1618  }
1619 
1620  /**
1621  * STEP 6: Now, all ingredients are ready to set up the actual
1622  * local to global mapping.
1623  *
1624  * The remainder of the map consists of the element-interior
1625  * degrees of freedom. This leads to the block-diagonal submatrix
1626  * as each element-interior mode is globally orthogonal to modes
1627  * in all other elements.
1628  */
1629  cnt = 0;
1630 
1631  // Loop over all the elements in the domain
1632  for(i = 0; i < locExpVector.size(); ++i)
1633  {
1634  exp = locExpVector[i];
1635  cnt = locExp.GetCoeff_Offset(i);
1636  for(j = 0; j < exp->GetNverts(); ++j)
1637  {
1638  meshVertId = exp->GetGeom()->GetVid(j);
1639 
1640  // Set the global DOF for vertex j of element i
1641  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1642  graphVertOffset[graph[0][meshVertId]];
1643  }
1644 
1645  for(j = 0; j < exp->GetNedges(); ++j)
1646  {
1647  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1648  edgeOrient = exp->GetGeom()->GetEorient(j);
1649  meshEdgeId = exp->GetGeom()->GetEid(j);
1650 
1651  pIt = periodicEdges.find(meshEdgeId);
1652 
1653  // See if this edge is periodic. If it is, then we map all
1654  // edges to the one with lowest ID, and align all
1655  // coefficients to this edge orientation.
1656  if (pIt != periodicEdges.end())
1657  {
1658  pair<int, StdRegions::Orientation> idOrient =
1660  meshEdgeId, edgeOrient, pIt->second);
1661  edgeOrient = idOrient.second;
1662  }
1663 
1664  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1665 
1666  // Set the global DOF's for the interior modes of edge j
1667  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1668  {
1669  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1670  graphVertOffset[graph[1][meshEdgeId]]+k;
1671  }
1672  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1673  {
1674  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1675  graphVertOffset[graph[1][meshEdgeId]];
1676  }
1677 
1678  // Fill the sign vector if required
1679  if(m_signChange)
1680  {
1681  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1682  {
1683  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1684  }
1685  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1686  {
1687  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = 0.0;
1688  }
1689  }
1690  }
1691 
1692  for(j = 0; j < exp->GetNfaces(); ++j)
1693  {
1694  faceOrient = exp->GetGeom()->GetForient(j);
1695  meshFaceId = exp->GetGeom()->GetFid(j);
1696 
1697  pIt = periodicFaces.find(meshFaceId);
1698 
1699  if (pIt != periodicFaces.end() &&
1700  meshFaceId == min(meshFaceId, pIt->second[0].id))
1701  {
1702  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1703  }
1704 
1705  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1706 
1707  // Set the global DOF's for the interior modes of face j
1708  exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
1709  switch(faceType[meshFaceId])
1710  {
1712  {
1713  int kLoc=0;
1714  k = 0;
1715  for( q = 2; q < numModes1; q++)
1716  {
1717  for( p = 2; p < numModes0; p++)
1718  {
1719  if( (p < faceModes[0][meshFaceId]) &&
1720  (q < faceModes[1][meshFaceId]))
1721  {
1722  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1723  graphVertOffset[graph[2][meshFaceId]]+k;
1724  if(m_signChange)
1725  {
1726  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1727  (NekDouble) faceInteriorSign[kLoc];
1728  }
1729  k++;
1730  }
1731  else
1732  {
1733  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1734  graphVertOffset[graph[2][meshFaceId]];
1735  if(m_signChange)
1736  {
1737  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1738  }
1739  }
1740  kLoc++;
1741  }
1742  }
1743  }
1744  break;
1746  {
1747  int kLoc=0;
1748  k = 0;
1749  for( p = 2; p < numModes0; p++)
1750  {
1751  for( q = 1; q < numModes1-p; q++)
1752  {
1753  if( (p < faceModes[0][meshFaceId]) &&
1754  (p+q < faceModes[1][meshFaceId]))
1755  {
1756  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1757  graphVertOffset[graph[2][meshFaceId]]+k;
1758  if(m_signChange)
1759  {
1760  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1761  (NekDouble) faceInteriorSign[kLoc];
1762  }
1763  k++;
1764  }
1765  else
1766  {
1767  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1768  graphVertOffset[graph[2][meshFaceId]];
1769  if(m_signChange)
1770  {
1771  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1772  }
1773  }
1774  kLoc++;
1775  }
1776  }
1777  }
1778  break;
1779  default:
1780  ASSERTL0(false,"Shape not recognised");
1781  break;
1782  }
1783  }
1784  }
1785 
1786  // Set up the mapping for the boundary conditions
1787  cnt = 0;
1788  int offset = 0;
1789  for(i = 0; i < bndCondExp.num_elements(); i++)
1790  {
1791  set<int> foundExtraVerts, foundExtraEdges;
1792  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1793  {
1794  bndExp = bndCondExp[i]->GetExp(j);
1795  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1796  for(k = 0; k < bndExp->GetNverts(); k++)
1797  {
1798  meshVertId = bndExp->GetGeom()->GetVid(k);
1799  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1800 
1801  if (bndConditions[i]->GetBoundaryConditionType() !=
1803  {
1804  continue;
1805  }
1806 
1807  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1808  if (iter != extraDirVerts.end() &&
1809  foundExtraVerts.count(meshVertId) == 0)
1810  {
1811  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1812  bndExp->GetVertexMap(k);
1813  int gid = graphVertOffset[
1814  graph[0][meshVertId]];
1815  ExtraDirDof t(loc, gid, 1.0);
1816  m_extraDirDofs[i].push_back(t);
1817  foundExtraVerts.insert(meshVertId);
1818  }
1819  }
1820 
1821  for(k = 0; k < bndExp->GetNedges(); k++)
1822  {
1823  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1824  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1825  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1826 
1827  pIt = periodicEdges.find(meshEdgeId);
1828 
1829  // See if this edge is periodic. If it is, then we map
1830  // all edges to the one with lowest ID, and align all
1831  // coefficients to this edge orientation.
1832  if (pIt != periodicEdges.end())
1833  {
1834  pair<int, StdRegions::Orientation> idOrient =
1836  meshEdgeId, edgeOrient, pIt->second);
1837  edgeOrient = idOrient.second;
1838  }
1839 
1840  bndExp->GetEdgeInteriorMap(
1841  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1842 
1843  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1844  {
1845  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1846  graphVertOffset[graph[1][meshEdgeId]]+l;
1847  }
1848  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1849  {
1850  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1851  graphVertOffset[graph[1][meshEdgeId]];
1852  }
1853 
1854  // Fill the sign vector if required
1855  if(m_signChange)
1856  {
1857  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1858  {
1859  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1860  }
1861  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1862  {
1863  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = 0.0;
1864  }
1865  }
1866 
1867  if (bndConditions[i]->GetBoundaryConditionType() !=
1869  {
1870  continue;
1871  }
1872 
1873  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1874  if (iter != extraDirEdges.end() &&
1875  foundExtraEdges.count(meshEdgeId) == 0 &&
1876  nEdgeInteriorCoeffs > 0)
1877  {
1878  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1879  {
1880  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1881  edgeInteriorMap[l];
1882  int gid = graphVertOffset[
1883  graph[1][meshEdgeId]]+l;
1884  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1885  m_extraDirDofs[i].push_back(t);
1886  }
1887  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1888  {
1889  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1890  edgeInteriorMap[l];
1891  int gid = graphVertOffset[
1892  graph[1][meshEdgeId]];
1893  ExtraDirDof t(loc, gid, 0.0);
1894  m_extraDirDofs[i].push_back(t);
1895  }
1896  foundExtraEdges.insert(meshEdgeId);
1897  }
1898  }
1899 
1900  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1901  intDofCnt = 0;
1902  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1903  {
1904  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1905  {
1907  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1908  intDofCnt++;
1909  }
1910  }
1911  }
1912  offset += bndCondExp[i]->GetNcoeffs();
1913  }
1914 
1915  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1916  m_numGlobalBndCoeffs = globalId;
1917 
1918 
1919  /*
1920  * The boundary condition mapping is generated from the same vertex
1921  * renumbering.
1922  */
1923  cnt=0;
1924  for(i = 0; i < m_numLocalCoeffs; ++i)
1925  {
1926  if(m_localToGlobalMap[i] == -1)
1927  {
1928  m_localToGlobalMap[i] = globalId++;
1929  }
1930  else
1931  {
1932  if(m_signChange)
1933  {
1935  }
1937  }
1938  }
1939 
1940  m_numGlobalCoeffs = globalId;
1941 
1942  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1943 
1944  // Since we can now have multiple entries to m_extraDirDofs due to
1945  // periodic boundary conditions we make a call to work out the
1946  // multiplicity of all entries and invert (Need to be after
1947  // SetupUniversalC0ContMap)
1948  Array<OneD, NekDouble> valence(m_numGlobalBndCoeffs,0.0);
1949 
1950  // Fill in Dirichlet coefficients that are to be sent to other
1951  // processors with a value of 1
1952  map<int, vector<ExtraDirDof> >::iterator Tit;
1953 
1954  // Generate valence for extraDirDofs
1955  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1956  {
1957  for (i = 0; i < Tit->second.size(); ++i)
1958  {
1959  valence[Tit->second[i].get<1>()] = 1.0;
1960  }
1961  }
1962 
1963  UniversalAssembleBnd(valence);
1964 
1965  // Set third argument of tuple to inverse of valence.
1966  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1967  {
1968  for (i = 0; i < Tit->second.size(); ++i)
1969  {
1970  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1971  }
1972  }
1973 
1974  // Set up the local to global map for the next level when using
1975  // multi-level static condensation
1979  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1980  {
1981  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1982  {
1983  Array<OneD, int> vwgts_perm(
1984  graph[0].size() + graph[1].size() + graph[2].size()
1985  - firstNonDirGraphVertId);
1986 
1987  for (i = 0; i < locExpVector.size(); ++i)
1988  {
1989  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1990 
1991  for (j = 0; j < exp->GetNverts(); ++j)
1992  {
1993  meshVertId = exp->GetGeom()->GetVid(j);
1994 
1995  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1996  {
1997  vwgts_perm[graph[0][meshVertId] -
1998  firstNonDirGraphVertId] =
1999  dofs[0][meshVertId];
2000  }
2001  }
2002 
2003  for (j = 0; j < exp->GetNedges(); ++j)
2004  {
2005  meshEdgeId = exp->GetGeom()->GetEid(j);
2006 
2007  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
2008  {
2009  vwgts_perm[graph[1][meshEdgeId] -
2010  firstNonDirGraphVertId] =
2011  dofs[1][meshEdgeId];
2012  }
2013  }
2014 
2015  for (j = 0; j < exp->GetNfaces(); ++j)
2016  {
2017  meshFaceId = exp->GetGeom()->GetFid(j);
2018 
2019  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
2020  {
2021  vwgts_perm[graph[2][meshFaceId] -
2022  firstNonDirGraphVertId] =
2023  dofs[2][meshFaceId];
2024  }
2025  }
2026  }
2027 
2028  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
2030  AllocateSharedPtr(this, bottomUpGraph);
2031  }
2032  }
2033 
2034  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
2035  m_localToGlobalMap.end());
2036 
2037  // Add up hash values if parallel
2038  int hash = m_hash;
2039  vComm->AllReduce(hash, LibUtilities::ReduceSum);
2040  m_hash = hash;
2041 
2044  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:120
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:316
Principle Modified Functions .
Definition: BasisType.h:51
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:308
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:239
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:222
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:779
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
Principle Modified Functions .
Definition: BasisType.h:49
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
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:397
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:311
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:390
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:365
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:356
boost::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:54
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
Principle Modified Functions .
Definition: BasisType.h:50
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
double NekDouble
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:392
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:350
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:139
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:354
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector< std::map< int, int > > DofGraph
Definition: AssemblyMapCG.h:56
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:314
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:386
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:352
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:79
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:305
void UniversalAssembleBnd(Array< OneD, NekDouble > &pGlobal) const
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
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...
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:388
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
Nektar::MultiRegions::AssemblyMapCG::~AssemblyMapCG ( )
virtual

Destructor.

Definition at line 2049 of file AssemblyMapCG.cpp.

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

2050  {
2053  }
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:222

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

2528  {
2529  int i,j;
2530  int cnt = 0;
2531  int locSize;
2532  int maxId;
2533  int minId;
2534  int bwidth = -1;
2535  for(i = 0; i < m_numPatches; ++i)
2536  {
2538  maxId = -1;
2539  minId = m_numLocalCoeffs+1;
2540  for(j = 0; j < locSize; j++)
2541  {
2543  {
2544  if(m_localToGlobalMap[cnt+j] > maxId)
2545  {
2546  maxId = m_localToGlobalMap[cnt+j];
2547  }
2548 
2549  if(m_localToGlobalMap[cnt+j] < minId)
2550  {
2551  minId = m_localToGlobalMap[cnt+j];
2552  }
2553  }
2554  }
2555  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2556 
2557  cnt+=locSize;
2558  }
2559 
2560  m_fullSystemBandWidth = bwidth;
2561  }
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
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:390
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:392
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:388
int Nektar::MultiRegions::AssemblyMapCG::CreateGraph ( const ExpList locExp,
const BndCondExp bndCondExp,
const Array< OneD, const BndCond > &  bndConditions,
const bool  checkIfSystemSingular,
const PeriodicMap periodicVerts,
const PeriodicMap periodicEdges,
const PeriodicMap periodicFaces,
DofGraph graph,
BottomUpSubStructuredGraphSharedPtr bottomUpGraph,
std::set< int > &  extraDirVerts,
std::set< int > &  extraDirEdges,
int &  firstNonDirGraphVertId,
int &  nExtraDirichlet,
int  mdswitch = 1 
)
protected
  • Count verts, edges, face and add up edges and face sizes
  • Periodic vertices
  • Periodic edges
  • Periodic faces

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

Definition at line 81 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(), CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceMin, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

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

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

Definition at line 108 of file AssemblyMapCG.h.

References m_extraDirDofs.

109  {
110  return m_extraDirDofs;
111  }
std::map< int, std::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 2165 of file AssemblyMapCG.cpp.

References Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetSession(), 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, m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Gs::Unique(), and Vmath::Zero().

Referenced by AssemblyMapCG().

2170  {
2172  int nVert = 0;
2173  int nEdge = 0;
2174  int nFace = 0;
2175  int maxEdgeDof = 0;
2176  int maxFaceDof = 0;
2177  int maxIntDof = 0;
2178  int dof = 0;
2179  int cnt;
2180  int i,j,k,l;
2181  int meshVertId;
2182  int meshEdgeId;
2183  int meshFaceId;
2184  int elementId;
2185  int vGlobalId;
2186  int maxBndGlobalId = 0;
2187  StdRegions::Orientation edgeOrient;
2188  StdRegions::Orientation faceOrient;
2189  Array<OneD, unsigned int> edgeInteriorMap;
2190  Array<OneD, int> edgeInteriorSign;
2191  Array<OneD, unsigned int> faceInteriorMap;
2192  Array<OneD, int> faceInteriorSign;
2193  Array<OneD, unsigned int> interiorMap;
2194  PeriodicMap::const_iterator pIt;
2195 
2196  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
2197  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
2198  const bool verbose = locExp.GetSession()->DefinesCmdLineArgument("verbose");
2199 
2204 
2205  // Loop over all the elements in the domain to gather mesh data
2206  for(i = 0; i < locExpVector.size(); ++i)
2207  {
2208  exp = locExpVector[i];
2209  nVert += exp->GetNverts();
2210  nEdge += exp->GetNedges();
2211  nFace += exp->GetNfaces();
2212  // Loop over all edges (and vertices) of element i
2213  for(j = 0; j < exp->GetNedges(); ++j)
2214  {
2215  dof = exp->GetEdgeNcoeffs(j)-2;
2216  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2217  }
2218  for(j = 0; j < exp->GetNfaces(); ++j)
2219  {
2220  dof = exp->GetFaceIntNcoeffs(j);
2221  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2222  }
2223  exp->GetInteriorMap(interiorMap);
2224  dof = interiorMap.num_elements();
2225  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2226  }
2227 
2228  // Tell other processes about how many dof we have
2229  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2230  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2231  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2232  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2233  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2234  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2235 
2236  // Assemble global to universal mapping for this process
2237  for(i = 0; i < locExpVector.size(); ++i)
2238  {
2239  exp = locExpVector[i];
2240  cnt = locExp.GetCoeff_Offset(i);
2241 
2242  // Loop over all vertices of element i
2243  for(j = 0; j < exp->GetNverts(); ++j)
2244  {
2245  meshVertId = exp->GetGeom()->GetVid(j);
2246  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2247 
2248  pIt = perVerts.find(meshVertId);
2249  if (pIt != perVerts.end())
2250  {
2251  for (k = 0; k < pIt->second.size(); ++k)
2252  {
2253  meshVertId = min(meshVertId, pIt->second[k].id);
2254  }
2255  }
2256 
2257  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2258  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2259  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2260  }
2261 
2262  // Loop over all edges of element i
2263  for(j = 0; j < exp->GetNedges(); ++j)
2264  {
2265  meshEdgeId = exp->GetGeom()->GetEid(j);
2266  pIt = perEdges.find(meshEdgeId);
2267  edgeOrient = exp->GetGeom()->GetEorient(j);
2268 
2269  if (pIt != perEdges.end())
2270  {
2271  pair<int, StdRegions::Orientation> idOrient =
2273  meshEdgeId, edgeOrient, pIt->second);
2274  meshEdgeId = idOrient.first;
2275  edgeOrient = idOrient.second;
2276  }
2277 
2278  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
2279  dof = exp->GetEdgeNcoeffs(j)-2;
2280 
2281  // Set the global DOF's for the interior modes of edge j
2282  // for varP, ignore modes with sign == 0
2283  for(k = 0, l = 0; k < dof; ++k)
2284  {
2285  if (m_signChange)
2286  {
2287  if (m_localToGlobalSign[cnt+edgeInteriorMap[k]]==0)
2288  {
2289  continue;
2290  }
2291  }
2292  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2293  m_globalToUniversalMap[vGlobalId]
2294  = nVert + meshEdgeId * maxEdgeDof + l + 1;
2295  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2296  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2297  l++;
2298  }
2299  }
2300 
2301  // Loop over all faces of element i
2302  for(j = 0; j < exp->GetNfaces(); ++j)
2303  {
2304  faceOrient = exp->GetGeom()->GetForient(j);
2305 
2306  meshFaceId = exp->GetGeom()->GetFid(j);
2307 
2308  pIt = perFaces.find(meshFaceId);
2309  if (pIt != perFaces.end())
2310  {
2311  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2312  {
2313  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
2314  }
2315  meshFaceId = min(meshFaceId, pIt->second[0].id);
2316  }
2317 
2318 
2319  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2320  dof = exp->GetFaceIntNcoeffs(j);
2321 
2322  for(k = 0, l = 0; k < dof; ++k)
2323  {
2324  if (m_signChange)
2325  {
2326  if (m_localToGlobalSign[cnt+faceInteriorMap[k]]==0)
2327  {
2328  continue;
2329  }
2330  }
2331  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2332  m_globalToUniversalMap[vGlobalId]
2333  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2334  + l + 1;
2335  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2336 
2337  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2338  l++;
2339  }
2340  }
2341 
2342  // Add interior DOFs to complete universal numbering
2343  exp->GetInteriorMap(interiorMap);
2344  dof = interiorMap.num_elements();
2345  elementId = (exp->GetGeom())->GetGlobalID();
2346  for (k = 0; k < dof; ++k)
2347  {
2348  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2349  m_globalToUniversalMap[vGlobalId]
2350  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2351  }
2352  }
2353 
2354  // Set up the GSLib universal assemble mapping
2355  // Internal DOF do not participate in any data
2356  // exchange, so we keep these set to the special GSLib id=0 so
2357  // they are ignored.
2358  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2359  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2360  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2361  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2362  {
2363  tmp[i] = m_globalToUniversalMap[i];
2364  }
2365 
2366  m_gsh = Gs::Init(tmp, vCommRow, verbose);
2367  m_bndGsh = Gs::Init(tmp2, vCommRow, verbose);
2368  Gs::Unique(tmp, vCommRow);
2369  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2370  {
2371  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2372  }
2373  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2374  {
2375  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2376  }
2377  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:316
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:308
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
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.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
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:200
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:360
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:362
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
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...
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
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 2687 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().

2690  {
2691  Array<OneD, const NekDouble> local;
2692  if(global.data() == loc.data())
2693  {
2694  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2695  }
2696  else
2697  {
2698  local = loc; // create reference
2699  }
2700  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2701 
2702  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2703 
2704  if(m_signChange)
2705  {
2706  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2707  }
2708  else
2709  {
2710  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2711  }
2712  UniversalAssemble(global);
2713  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
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:709
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
void Nektar::MultiRegions::AssemblyMapCG::v_Assemble ( const NekVector< NekDouble > &  loc,
NekVector< NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2715 of file AssemblyMapCG.cpp.

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

2718  {
2719  Assemble(loc.GetPtr(),global.GetPtr());
2720  }
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 2784 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2785  {
2786  return m_extraDirEdges;
2787  }
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 2744 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2745  {
2746  return m_fullSystemBandWidth;
2747  }
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 2569 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2570  {
2571  return m_globalToUniversalMap[i];
2572  }
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 2586 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2587  {
2588  return m_globalToUniversalMap;
2589  }
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 2574 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2575  {
2576  return m_globalToUniversalMapUnique[i];
2577  }
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 2592 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2593  {
2595  }
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 2564 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2565  {
2566  return m_localToGlobalMap[i];
2567  }
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 2580 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2581  {
2582  return m_localToGlobalMap;
2583  }
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 2597 of file AssemblyMapCG.cpp.

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

2599  {
2600  if(m_signChange)
2601  {
2602  return m_localToGlobalSign[i];
2603  }
2604  else
2605  {
2606  return 1.0;
2607  }
2608  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
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 2610 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2611  {
2612  return m_localToGlobalSign;
2613  }
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 2764 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2765  {
2766  return m_numDirEdges;
2767  }
int m_numDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2769 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2770  {
2771  return m_numDirFaces;
2772  }
int m_numDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2754 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2755  {
2756  return m_numNonDirEdgeModes;
2757  }
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 2774 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2775  {
2776  return m_numNonDirEdges;
2777  }
int m_numNonDirEdges
Number of Dirichlet edges.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2759 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2760  {
2761  return m_numNonDirFaceModes;
2762  }
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 2779 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2780  {
2781  return m_numNonDirFaces;
2782  }
int m_numNonDirFaces
Number of Dirichlet faces.
int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2749 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2750  {
2751  return m_numNonDirVertexModes;
2752  }
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 2655 of file AssemblyMapCG.cpp.

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

2658  {
2659  Array<OneD, const NekDouble> glo;
2660  if(global.data() == loc.data())
2661  {
2662  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2663  }
2664  else
2665  {
2666  glo = global; // create reference
2667  }
2668 
2669 
2670  if(m_signChange)
2671  {
2672  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2673  }
2674  else
2675  {
2676  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2677  }
2678  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:644
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
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 2680 of file AssemblyMapCG.cpp.

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

2683  {
2684  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2685  }
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 2387 of file AssemblyMapCG.cpp.

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

2389  {
2390  AssemblyMapCGSharedPtr returnval;
2391 
2392  int i, j;
2393  int nverts = 0;
2394  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2395  = locexp.GetExp();
2396  int nelmts = exp->size();
2397  const bool verbose = locexp.GetSession()->DefinesCmdLineArgument("verbose");
2398 
2399  // Get Default Map and turn off any searched values.
2400  returnval = MemoryManager<AssemblyMapCG>
2402  returnval->m_solnType = solnType;
2403  returnval->m_preconType = eNull;
2404  returnval->m_maxStaticCondLevel = 0;
2405  returnval->m_signChange = false;
2406  returnval->m_comm = m_comm;
2407 
2408  // Count the number of vertices
2409  for (i = 0; i < nelmts; ++i)
2410  {
2411  nverts += (*exp)[i]->GetNverts();
2412  }
2413 
2414  returnval->m_numLocalCoeffs = nverts;
2415  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2416 
2417  // Store original global ids in this map
2418  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2419 
2420  int cnt = 0;
2421  int cnt1 = 0;
2422  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2423 
2424  // Set up local to global map;
2425  for (i = 0; i < nelmts; ++i)
2426  {
2427  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2428  {
2429  returnval->m_localToGlobalMap[cnt] =
2430  returnval->m_localToGlobalBndMap[cnt] =
2431  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2432  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2433 
2434  // Set up numLocalDirBndCoeffs
2435  if ((returnval->m_localToGlobalMap[cnt]) <
2437  {
2438  returnval->m_numLocalDirBndCoeffs++;
2439  }
2440  cnt++;
2441  }
2442  cnt1 += (*exp)[i]->GetNcoeffs();
2443  }
2444 
2445  cnt = 0;
2446  // Reset global numbering and count number of dofs
2447  for (i = 0; i < m_numGlobalCoeffs; ++i)
2448  {
2449  if (GlobCoeffs[i] != -1)
2450  {
2451  GlobCoeffs[i] = cnt++;
2452  }
2453  }
2454 
2455  // Set up number of globalCoeffs;
2456  returnval->m_numGlobalCoeffs = cnt;
2457 
2458  // Set up number of global Dirichlet boundary coefficients
2459  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2460  {
2461  if (GlobCoeffs[i] != -1)
2462  {
2463  returnval->m_numGlobalDirBndCoeffs++;
2464  }
2465  }
2466 
2467  // Set up global to universal map
2468  if (m_globalToUniversalMap.num_elements())
2469  {
2471  = m_session->GetComm()->GetRowComm();
2472  int nglocoeffs = returnval->m_numGlobalCoeffs;
2473  returnval->m_globalToUniversalMap
2474  = Array<OneD, int> (nglocoeffs);
2475  returnval->m_globalToUniversalMapUnique
2476  = Array<OneD, int> (nglocoeffs);
2477 
2478  // Reset local to global map and setup universal map
2479  for (i = 0; i < nverts; ++i)
2480  {
2481  cnt = returnval->m_localToGlobalMap[i];
2482  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2483 
2484  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2486  }
2487 
2488  Nektar::Array<OneD, long> tmp(nglocoeffs);
2489  Vmath::Zero(nglocoeffs, tmp, 1);
2490  for (unsigned int i = 0; i < nglocoeffs; ++i)
2491  {
2492  tmp[i] = returnval->m_globalToUniversalMap[i];
2493  }
2494  returnval->m_gsh = Gs::Init(tmp, vCommRow, verbose);
2495  Gs::Unique(tmp, vCommRow);
2496  for (unsigned int i = 0; i < nglocoeffs; ++i)
2497  {
2498  returnval->m_globalToUniversalMapUnique[i]
2499  = (tmp[i] >= 0 ? 1 : 0);
2500  }
2501  }
2502  else // not sure this option is ever needed.
2503  {
2504  for (i = 0; i < nverts; ++i)
2505  {
2506  cnt = returnval->m_localToGlobalMap[i];
2507  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2508  }
2509  }
2510 
2511  return returnval;
2512  }
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:308
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
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:200
No Solution type specified.
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:305
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
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,
bool  useComm 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2615 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().

2619  {
2620  Array<OneD, const NekDouble> local;
2621  if(global.data() == loc.data())
2622  {
2623  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2624  }
2625  else
2626  {
2627  local = loc; // create reference
2628  }
2629 
2630 
2631  if(m_signChange)
2632  {
2633  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2634  }
2635  else
2636  {
2637  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2638  }
2639 
2640  // ensure all values are unique by calling a max
2641  if(useComm)
2642  {
2643  Gs::Gather(global, Gs::gs_max, m_gsh);
2644  }
2645  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
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:239
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
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:673
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,
bool  useComm 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2647 of file AssemblyMapCG.cpp.

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

2651  {
2652  LocalToGlobal(loc.GetPtr(),global.GetPtr(),useComm);
2653  }
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) 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 2722 of file AssemblyMapCG.cpp.

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

2724  {
2725  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2726  }
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:239
void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( NekVector< NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2728 of file AssemblyMapCG.cpp.

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

2730  {
2731  UniversalAssemble(pGlobal.GetPtr());
2732  }
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 2734 of file AssemblyMapCG.cpp.

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

2737  {
2738  Array<OneD, NekDouble> tmp(offset);
2739  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2740  UniversalAssemble(pGlobal);
2741  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2742  }
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:1061

Member Data Documentation

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

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

Definition at line 148 of file AssemblyMapCG.h.

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

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

Extra dirichlet edges in parallel.

Definition at line 141 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetExtraDirEdges().

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

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

Definition at line 119 of file AssemblyMapCG.h.

Referenced by CalculateFullSystemBandWidth(), and v_GetFullSystemBandWidth().

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

Integer map of process coeffs to universal space.

Definition at line 121 of file AssemblyMapCG.h.

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

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

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

Definition at line 123 of file AssemblyMapCG.h.

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

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