Nektar++
Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | List of all members
Nektar::MultiRegions::AssemblyMapCG Class Reference

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

#include <AssemblyMapCG.h>

Inheritance diagram for Nektar::MultiRegions::AssemblyMapCG:
[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 ()
 
std::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 ExpListSharedPtrBndCondExp
 
typedef Array< OneD, const SpatialDomains::BoundaryConditionShPtrBndCond
 

Detailed Description

Constructs mappings for the C0 scalar continuous Galerkin formulation.

Mappings are created for three possible global solution types:

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

Definition at line 70 of file AssemblyMapCG.h.

Member Typedef Documentation

◆ BndCond

Definition at line 74 of file AssemblyMapCG.h.

◆ BndCondExp

Definition at line 72 of file AssemblyMapCG.h.

Constructor & Destructor Documentation

◆ AssemblyMapCG() [1/2]

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

Default constructor.

Definition at line 71 of file AssemblyMapCG.cpp.

References m_maxStaticCondLevel.

73  :
74  AssemblyMap(pSession,variable)
75  {
76  pSession->LoadParameter(
77  "MaxStaticCondLevel", m_maxStaticCondLevel, 100);
78  }
int m_maxStaticCondLevel
Maximum static condensation level.
AssemblyMap()
Default constructor.
Definition: AssemblyMap.cpp:80

◆ AssemblyMapCG() [2/2]

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 1302 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::LibUtilities::eModifiedPyr_C, Nektar::SpatialDomains::ePeriodic, 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(), Gs::gs_min, Nektar::hash_range(), Gs::Init(), CG_Iterations::loc, 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().

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

◆ ~AssemblyMapCG()

Nektar::MultiRegions::AssemblyMapCG::~AssemblyMapCG ( )
virtual

Destructor.

Definition at line 2064 of file AssemblyMapCG.cpp.

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

2065  {
2068  }
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:226

Member Function Documentation

◆ CalculateFullSystemBandWidth()

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

2542  {
2543  int i,j;
2544  int cnt = 0;
2545  int locSize;
2546  int maxId;
2547  int minId;
2548  int bwidth = -1;
2549  for(i = 0; i < m_numPatches; ++i)
2550  {
2552  maxId = -1;
2553  minId = m_numLocalCoeffs+1;
2554  for(j = 0; j < locSize; j++)
2555  {
2557  {
2558  if(m_localToGlobalMap[cnt+j] > maxId)
2559  {
2560  maxId = m_localToGlobalMap[cnt+j];
2561  }
2562 
2563  if(m_localToGlobalMap[cnt+j] < minId)
2564  {
2565  minId = m_localToGlobalMap[cnt+j];
2566  }
2567  }
2568  }
2569  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2570 
2571  cnt+=locSize;
2572  }
2573 
2574  m_fullSystemBandWidth = bwidth;
2575  }
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:332
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:389
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:319
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:391
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:387

◆ CreateGraph()

int Nektar::MultiRegions::AssemblyMapCG::CreateGraph ( const ExpList locExp,
const BndCondExp bndCondExp,
const Array< OneD, const BndCond > &  bndConditions,
const bool  checkIfSystemSingular,
const PeriodicMap periodicVerts,
const PeriodicMap periodicEdges,
const PeriodicMap periodicFaces,
DofGraph graph,
BottomUpSubStructuredGraphSharedPtr bottomUpGraph,
std::set< int > &  extraDirVerts,
std::set< int > &  extraDirEdges,
int &  firstNonDirGraphVertId,
int &  nExtraDirichlet,
int  mdswitch = 1 
)
protected
  • Count verts, edges, face and add up edges and face sizes
  • Periodic vertices
  • Periodic edges
  • Periodic faces

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

Definition at line 80 of file AssemblyMapCG.cpp.

References ASSERTL0, ASSERTL1, Nektar::MultiRegions::CuthillMckeeReordering(), Nektar::MultiRegions::eDirectFullMatrix, Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::MultiRegions::eDirectStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::MultiRegions::eIterativeFull, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::MultiRegions::eIterativeStaticCond, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, 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::GlobalSysSolnTypeMap, Gs::gs_add, Vmath::Imax(), Gs::Init(), 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().

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

◆ GetExtraDirDofs()

std::map<int, std::vector<ExtraDirDof> >& Nektar::MultiRegions::AssemblyMapCG::GetExtraDirDofs ( )
inline

Definition at line 107 of file AssemblyMapCG.h.

References m_extraDirDofs.

108  {
109  return m_extraDirDofs;
110  }
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...

◆ SetUpUniversalC0ContMap()

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

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

◆ v_Assemble() [1/2]

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

References Vmath::Assmb(), CG_Iterations::loc, 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().

2704  {
2705  Array<OneD, const NekDouble> local;
2706  if(global.data() == loc.data())
2707  {
2708  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2709  }
2710  else
2711  {
2712  local = loc; // create reference
2713  }
2714  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2715 
2716  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2717 
2718  if(m_signChange)
2719  {
2720  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2721  }
2722  else
2723  {
2724  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2725  }
2726  UniversalAssemble(global);
2727  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:346
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:332
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero&#39;d first.
Definition: Vmath.cpp:712
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:343
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376

◆ v_Assemble() [2/2]

void Nektar::MultiRegions::AssemblyMapCG::v_Assemble ( const NekVector< NekDouble > &  loc,
NekVector< NekDouble > &  global 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2729 of file AssemblyMapCG.cpp.

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

2732  {
2733  Assemble(loc.GetPtr(),global.GetPtr());
2734  }
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:227

◆ v_GetExtraDirEdges()

const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetExtraDirEdges ( )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2798 of file AssemblyMapCG.cpp.

References m_extraDirEdges.

2799  {
2800  return m_extraDirEdges;
2801  }
Array< OneD, int > m_extraDirEdges
Extra dirichlet edges in parallel.

◆ v_GetFullSystemBandWidth()

int Nektar::MultiRegions::AssemblyMapCG::v_GetFullSystemBandWidth ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2758 of file AssemblyMapCG.cpp.

References m_fullSystemBandWidth.

2759  {
2760  return m_fullSystemBandWidth;
2761  }
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).

◆ v_GetGlobalToUniversalMap() [1/2]

int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2583 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2584  {
2585  return m_globalToUniversalMap[i];
2586  }
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.

◆ v_GetGlobalToUniversalMap() [2/2]

const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2600 of file AssemblyMapCG.cpp.

References m_globalToUniversalMap.

2601  {
2602  return m_globalToUniversalMap;
2603  }
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.

◆ v_GetGlobalToUniversalMapUnique() [1/2]

int Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2588 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2589  {
2590  return m_globalToUniversalMapUnique[i];
2591  }
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)

◆ v_GetGlobalToUniversalMapUnique() [2/2]

const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetGlobalToUniversalMapUnique ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2606 of file AssemblyMapCG.cpp.

References m_globalToUniversalMapUnique.

2607  {
2609  }
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)

◆ v_GetLocalToGlobalMap() [1/2]

int Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2578 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2579  {
2580  return m_localToGlobalMap[i];
2581  }
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.

◆ v_GetLocalToGlobalMap() [2/2]

const Array< OneD, const int > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalMap ( void  )
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2594 of file AssemblyMapCG.cpp.

References m_localToGlobalMap.

2595  {
2596  return m_localToGlobalMap;
2597  }
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.

◆ v_GetLocalToGlobalSign() [1/2]

NekDouble Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( const int  i) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2611 of file AssemblyMapCG.cpp.

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

2613  {
2614  if(m_signChange)
2615  {
2616  return m_localToGlobalSign[i];
2617  }
2618  else
2619  {
2620  return 1.0;
2621  }
2622  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:346
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.

◆ v_GetLocalToGlobalSign() [2/2]

const Array< OneD, NekDouble > & Nektar::MultiRegions::AssemblyMapCG::v_GetLocalToGlobalSign ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2624 of file AssemblyMapCG.cpp.

References m_localToGlobalSign.

2625  {
2626  return m_localToGlobalSign;
2627  }
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.

◆ v_GetNumDirEdges()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2778 of file AssemblyMapCG.cpp.

References m_numDirEdges.

2779  {
2780  return m_numDirEdges;
2781  }
int m_numDirEdges
Number of Dirichlet edges.

◆ v_GetNumDirFaces()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2783 of file AssemblyMapCG.cpp.

References m_numDirFaces.

2784  {
2785  return m_numDirFaces;
2786  }
int m_numDirFaces
Number of Dirichlet faces.

◆ v_GetNumNonDirEdgeModes()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdgeModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2768 of file AssemblyMapCG.cpp.

References m_numNonDirEdgeModes.

2769  {
2770  return m_numNonDirEdgeModes;
2771  }
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.

◆ v_GetNumNonDirEdges()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirEdges ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2788 of file AssemblyMapCG.cpp.

References m_numNonDirEdges.

2789  {
2790  return m_numNonDirEdges;
2791  }
int m_numNonDirEdges
Number of Dirichlet edges.

◆ v_GetNumNonDirFaceModes()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaceModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2773 of file AssemblyMapCG.cpp.

References m_numNonDirFaceModes.

2774  {
2775  return m_numNonDirFaceModes;
2776  }
int m_numNonDirFaceModes
Number of non Dirichlet face modes.

◆ v_GetNumNonDirFaces()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirFaces ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2793 of file AssemblyMapCG.cpp.

References m_numNonDirFaces.

2794  {
2795  return m_numNonDirFaces;
2796  }
int m_numNonDirFaces
Number of Dirichlet faces.

◆ v_GetNumNonDirVertexModes()

int Nektar::MultiRegions::AssemblyMapCG::v_GetNumNonDirVertexModes ( ) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2763 of file AssemblyMapCG.cpp.

References m_numNonDirVertexModes.

2764  {
2765  return m_numNonDirVertexModes;
2766  }
int m_numNonDirVertexModes
Number of non Dirichlet vertex modes.

◆ v_GlobalToLocal() [1/2]

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

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

2672  {
2673  Array<OneD, const NekDouble> glo;
2674  if(global.data() == loc.data())
2675  {
2676  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2677  }
2678  else
2679  {
2680  glo = global; // create reference
2681  }
2682 
2683 
2684  if(m_signChange)
2685  {
2686  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2687  }
2688  else
2689  {
2690  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2691  }
2692  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:346
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:647
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:332
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.

◆ v_GlobalToLocal() [2/2]

void Nektar::MultiRegions::AssemblyMapCG::v_GlobalToLocal ( const NekVector< NekDouble > &  global,
NekVector< NekDouble > &  loc 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2694 of file AssemblyMapCG.cpp.

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

2697  {
2698  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2699  }
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:227
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const

◆ v_LinearSpaceMap()

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

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

◆ v_LocalToGlobal() [1/2]

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

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

2633  {
2634  Array<OneD, const NekDouble> local;
2635  if(global.data() == loc.data())
2636  {
2637  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2638  }
2639  else
2640  {
2641  local = loc; // create reference
2642  }
2643 
2644 
2645  if(m_signChange)
2646  {
2647  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2648  }
2649  else
2650  {
2651  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2652  }
2653 
2654  // ensure all values are unique by calling a max
2655  if(useComm)
2656  {
2657  Gs::Gather(global, Gs::gs_max, m_gsh);
2658  }
2659  }
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:346
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:245
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:332
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:676
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.

◆ v_LocalToGlobal() [2/2]

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

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

2665  {
2666  LocalToGlobal(loc.GetPtr(),global.GetPtr(),useComm);
2667  }
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:227

◆ v_UniversalAssemble() [1/3]

void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2736 of file AssemblyMapCG.cpp.

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

2738  {
2739  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2740  }
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:245

◆ v_UniversalAssemble() [2/3]

void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( NekVector< NekDouble > &  pGlobal) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2742 of file AssemblyMapCG.cpp.

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

2744  {
2745  UniversalAssemble(pGlobal.GetPtr());
2746  }
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:227

◆ v_UniversalAssemble() [3/3]

void Nektar::MultiRegions::AssemblyMapCG::v_UniversalAssemble ( Array< OneD, NekDouble > &  pGlobal,
int  offset 
) const
protectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2748 of file AssemblyMapCG.cpp.

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

2751  {
2752  Array<OneD, NekDouble> tmp(offset);
2753  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2754  UniversalAssemble(pGlobal);
2755  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2756  }
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:1064

Member Data Documentation

◆ m_extraDirDofs

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

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

◆ m_extraDirEdges

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

Extra dirichlet edges in parallel.

Definition at line 140 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetExtraDirEdges().

◆ m_fullSystemBandWidth

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

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

Definition at line 118 of file AssemblyMapCG.h.

Referenced by CalculateFullSystemBandWidth(), and v_GetFullSystemBandWidth().

◆ m_globalToUniversalMap

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

Integer map of process coeffs to universal space.

Definition at line 120 of file AssemblyMapCG.h.

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

◆ m_globalToUniversalMapUnique

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

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

Definition at line 122 of file AssemblyMapCG.h.

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

◆ m_localToGlobalMap

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

◆ m_localToGlobalSign

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

◆ m_maxStaticCondLevel

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

Maximum static condensation level.

Definition at line 144 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG().

◆ m_numDirEdges

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

Number of Dirichlet edges.

Definition at line 130 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirEdges().

◆ m_numDirFaces

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

Number of Dirichlet faces.

Definition at line 132 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirFaces().

◆ m_numLocalBndCondCoeffs

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

Number of local boundary condition coefficients.

Definition at line 138 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and CreateGraph().

◆ m_numLocDirBndCondDofs

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

Number of local boundary condition degrees of freedom.

Definition at line 142 of file AssemblyMapCG.h.

◆ m_numNonDirEdgeModes

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

Number of non Dirichlet edge modes.

Definition at line 126 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().

◆ m_numNonDirEdges

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

Number of Dirichlet edges.

Definition at line 134 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdges().

◆ m_numNonDirFaceModes

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

Number of non Dirichlet face modes.

Definition at line 128 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().

◆ m_numNonDirFaces

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

Number of Dirichlet faces.

Definition at line 136 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaces().

◆ m_numNonDirVertexModes

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

Number of non Dirichlet vertex modes.

Definition at line 124 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().