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::set< ExtraDirDof > & GetCopyLocalDirDofs ()
 
std::set< int > & GetParallelDirBndSign ()
 
- 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
 
void UniversalAbsMaxBnd (Array< OneD, NekDouble > &bndvals)
 
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 NekDouble > GetLocalToGlobalBndSign () const
 Retrieve the sign change for all local boundary modes. More...
 
const Array< OneD, const int > & GetBndCondCoeffsToLocalCoeffsMap ()
 Retrieves the local indices corresponding to the boundary expansion modes. More...
 
const Array< OneD, NekDouble > & GetBndCondCoeffsToLocalCoeffsSign ()
 Returns the modal sign associated with a given boundary expansion mode. More...
 
const Array< OneD, const int > & GetBndCondCoeffsToLocalTraceMap ()
 Retrieves the local indices corresponding to the boundary expansion modes to global trace. More...
 
int GetBndCondIDToGlobalTraceID (const int i)
 Returns the global index of the boundary trace giving the index on the boundary expansion. More...
 
const Array< OneD, const int > & GetBndCondIDToGlobalTraceID ()
 
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 Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset, bool UseComm=true) const
 
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool UseComm=true) const
 
void LocalToLocalBnd (const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locbnd) const
 
void LocalToLocalInt (const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locint) const
 
void LocalBndToLocal (const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) const
 
void LocalIntToLocal (const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) 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< AssemblyMap > LinearSpaceMap (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 PatchMapSharedPtr & GetPatchMapFromPrevLevel (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
 
std::string GetLinSysIterSolver () const
 
int GetLowestStaticCondLevel () const
 
void PatchLocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
 
void PatchGlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
 
void PatchAssemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) 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, NekDouble > m_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::set< ExtraDirDof > m_copyLocalDirDofs
 Set indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. More...
 
std::set< int > m_parallelDirBndSign
 Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change. 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 coeffs to global Boundary Dofs. More...
 
Array< OneD, NekDouble > m_localToGlobalBndSign
 Integer sign of local boundary coeffs to global space. More...
 
Array< OneD, int > m_localToLocalBndMap
 Integer map of local boundary coeffs to local boundary system numbering. More...
 
Array< OneD, int > m_localToLocalIntMap
 Integer map of local boundary coeffs to local interior system numbering. More...
 
Array< OneD, int > m_bndCondCoeffsToLocalCoeffsMap
 Integer map of bnd cond coeffs to local coefficients. More...
 
Array< OneD, NekDouble > m_bndCondCoeffsToLocalCoeffsSign
 Integer map of sign of bnd cond coeffs to local coefficients. More...
 
Array< OneD, int > m_bndCondCoeffsToLocalTraceMap
 Integer map of bnd cond coeff to local trace coeff. More...
 
Array< OneD, int > m_bndCondIDToGlobalTraceID
 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...
 
std::string m_linSysIterSolver
 Iterative solver: Conjugate Gradient, GMRES. More...
 
Gs::gs_data * m_gsh
 
Gs::gs_data * m_bndGsh
 
Gs::gs_data * m_dirBndGsh
 gs gather communication to impose Dirhichlet BCs. More...
 
int m_staticCondLevel
 The level of recursion in the case of multi-level static condensation. More...
 
int m_numPatches
 The number of patches (~elements) in the current level. More...
 
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
 The number of bnd dofs per patch. More...
 
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
 The number of int dofs per patch. More...
 
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
 Map from the patches of the previous level to the patches of the current level. More...
 
int m_lowestStaticCondLevel
 Lowest static condensation level. More...
 

Private Types

typedef Array< OneD, const ExpListSharedPtr > BndCondExp
 
typedef Array< OneD, const SpatialDomains::BoundaryConditionShPtr > BndCond
 

Detailed Description

Constructs mappings for the C0 scalar continuous Galerkin formulation.

Mappings are created for three possible global solution types:

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

Definition at line 70 of file AssemblyMapCG.h.

Member Typedef Documentation

◆ BndCond

typedef Array<OneD, const SpatialDomains::BoundaryConditionShPtr> Nektar::MultiRegions::AssemblyMapCG::BndCond
private

Definition at line 74 of file AssemblyMapCG.h.

◆ BndCondExp

typedef Array<OneD, const ExpListSharedPtr> Nektar::MultiRegions::AssemblyMapCG::BndCondExp
private

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.

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

References m_maxStaticCondLevel.

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

1343  : AssemblyMap(pSession, variable)
1344  {
1345  int i, j, k;
1346  int p, q, numModes0, numModes1;
1347  int cnt = 0;
1348  int meshVertId, meshEdgeId, meshEdgeId2, meshFaceId, meshFaceId2;
1349  int globalId;
1350  int nEdgeInteriorCoeffs;
1351  int firstNonDirGraphVertId;
1352  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1353  LocalRegions::ExpansionSharedPtr exp, bndExp;
1354  StdRegions::Orientation edgeOrient;
1355  StdRegions::Orientation faceOrient;
1356  Array<OneD, unsigned int> edgeInteriorMap;
1357  Array<OneD, int> edgeInteriorSign;
1358  Array<OneD, unsigned int> faceInteriorMap;
1359  Array<OneD, int> faceInteriorSign;
1360 
1361  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1362 
1363  bool verbose = m_session->DefinesCmdLineArgument("verbose");
1364 
1365  m_signChange = false;
1366 
1367  // Stores vertex, edge and face reordered vertices.
1368  DofGraph graph(3);
1369  DofGraph dofs(3);
1370  vector<map<int, int> > faceModes(2);
1371  map<int, LibUtilities::ShapeType> faceType;
1372 
1373  set<int> extraDirVerts, extraDirEdges;
1374  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
1375 
1376  // Construct list of number of degrees of freedom for each vertex,
1377  // edge and face.
1378  for (i = 0; i < locExpVector.size(); ++i)
1379  {
1380  exp = locExpVector[i];
1381 
1382  for(j = 0; j < exp->GetNverts(); ++j)
1383  {
1384  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1385  }
1386 
1387  for(j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1388  {
1389  int nEdgeInt;
1390  if(exp->GetGeom()->GetNumFaces())
1391  {
1392  nEdgeInt = exp->as<LocalRegions::Expansion3D>()->
1393  GetEdgeNcoeffs(j)-2;
1394  }
1395  else
1396  {
1397  nEdgeInt = exp->GetTraceNcoeffs(j)-2;
1398  }
1399 
1400  if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1401  {
1402  if (dofs[1][exp->GetGeom()->GetEid(j)] != nEdgeInt)
1403  {
1404  ASSERTL0( (exp->GetBasisType(0) == LibUtilities::eModified_A) ||
1405  (exp->GetBasisType(1) == LibUtilities::eModified_B) ||
1406  (exp->GetBasisType(2) == LibUtilities::eModified_C) ||
1407  (exp->GetBasisType(2) == LibUtilities::eModifiedPyr_C),
1408  "CG with variable order only available with modal expansion");
1409  }
1410  dofs[1][exp->GetGeom()->GetEid(j)] =
1411  min(dofs[1][exp->GetGeom()->GetEid(j)],
1412  nEdgeInt);
1413  }
1414  else
1415  {
1416  dofs[1][exp->GetGeom()->GetEid(j)] = nEdgeInt;
1417  }
1418  }
1419 
1420  for(j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1421  {
1422  faceOrient = exp->GetGeom()->GetForient(j);
1423  meshFaceId = exp->GetGeom()->GetFid(j);
1424  exp->GetTraceNumModes(j, numModes0, numModes1, faceOrient);
1425 
1426  if (faceModes[0].count(meshFaceId) > 0)
1427  {
1428  faceModes[0][meshFaceId] =
1429  min(faceModes[0][meshFaceId], numModes0);
1430 
1431  faceModes[1][meshFaceId] =
1432  min(faceModes[1][meshFaceId], numModes1);
1433  }
1434  else
1435  {
1436  faceModes[0][meshFaceId] = numModes0;
1437  faceModes[1][meshFaceId] = numModes1;
1438 
1439  // Get shape of this face
1440  SpatialDomains::Geometry3DSharedPtr geom;
1441  geom = std::dynamic_pointer_cast<SpatialDomains::
1442  Geometry3D> (exp->GetGeom());
1443  faceType[meshFaceId] =
1444  geom->GetFace(j)->GetShapeType();
1445  }
1446  }
1447  }
1448 
1449  // Add non-local periodic dofs to the map
1450  for (auto &pIt : periodicEdges)
1451  {
1452  for (i = 0; i < pIt.second.size(); ++i)
1453  {
1454  meshEdgeId2 = pIt.second[i].id;
1455  if (dofs[1].count(meshEdgeId2) == 0)
1456  {
1457  dofs[1][meshEdgeId2] = 1e6;
1458  }
1459  }
1460  }
1461  for (auto &pIt : periodicFaces)
1462  {
1463  for (i = 0; i < pIt.second.size(); ++i)
1464  {
1465  meshFaceId2 = pIt.second[i].id;
1466  if (faceModes[0].count(meshFaceId2) == 0)
1467  {
1468  faceModes[0][meshFaceId2] = 1e6;
1469  faceModes[1][meshFaceId2] = 1e6;
1470  }
1471  }
1472  }
1473 
1474  // Now use information from all partitions to determine the correct
1475  // size
1476 
1477  // edges
1478  Array<OneD, long> edgeId (dofs[1].size());
1479  Array<OneD, NekDouble> edgeDof (dofs[1].size());
1480  i = 0;
1481  for(auto &dofIt : dofs[1])
1482  {
1483  edgeId [i ] = dofIt.first + 1;
1484  edgeDof[i++] = (NekDouble) dofIt.second;
1485  }
1486  Gs::gs_data *tmp = Gs::Init(edgeId, vComm, verbose);
1487  Gs::Gather(edgeDof, Gs::gs_min, tmp);
1488  Gs::Finalise(tmp);
1489  for (i = 0; i < dofs[1].size(); i++)
1490  {
1491  dofs[1][edgeId[i]-1] = (int) (edgeDof[i]+0.5);
1492  }
1493  // Periodic edges
1494  for (auto &pIt : periodicEdges)
1495  {
1496  meshEdgeId = pIt.first;
1497  for (i = 0; i < pIt.second.size(); ++i)
1498  {
1499  meshEdgeId2 = pIt.second[i].id;
1500  if (dofs[1][meshEdgeId2] < dofs[1][meshEdgeId])
1501  {
1502  dofs[1][meshEdgeId] = dofs[1][meshEdgeId2];
1503  }
1504  }
1505  }
1506  // faces
1507  Array<OneD, long> faceId (faceModes[0].size());
1508  Array<OneD, NekDouble> faceP (faceModes[0].size());
1509  Array<OneD, NekDouble> faceQ (faceModes[0].size());
1510 
1511  i = 0;
1512  for(auto dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin();
1513  dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
1514  {
1515  faceId[i] = dofIt->first+1;
1516  faceP[i] = (NekDouble) dofIt->second;
1517  faceQ[i] = (NekDouble) dofIt2->second;
1518  }
1519  Gs::gs_data *tmp2 = Gs::Init(faceId, vComm, verbose);
1520  Gs::Gather(faceP, Gs::gs_min, tmp2);
1521  Gs::Gather(faceQ, Gs::gs_min, tmp2);
1522  Gs::Finalise(tmp2);
1523  for (i=0; i < faceModes[0].size(); i++)
1524  {
1525  faceModes[0][faceId[i]-1] = (int) (faceP[i]+0.5);
1526  faceModes[1][faceId[i]-1] = (int) (faceQ[i]+0.5);
1527  }
1528  // Periodic faces
1529  for (auto &pIt : periodicFaces)
1530  {
1531  meshFaceId = pIt.first;
1532  for (i = 0; i < pIt.second.size(); ++i)
1533  {
1534  meshFaceId2 = pIt.second[i].id;
1535  if (faceModes[0][meshFaceId2] < faceModes[0][meshFaceId])
1536  {
1537  faceModes[0][meshFaceId] = faceModes[0][meshFaceId2];
1538  }
1539  if (faceModes[1][meshFaceId2] < faceModes[1][meshFaceId])
1540  {
1541  faceModes[1][meshFaceId] = faceModes[1][meshFaceId2];
1542  }
1543  }
1544  }
1545  // Calculate number of dof in each face
1546  int P, Q;
1547  for (i=0; i < faceModes[0].size(); i++)
1548  {
1549  P = faceModes[0][faceId[i]-1];
1550  Q = faceModes[1][faceId[i]-1];
1551  if (faceType[faceId[i]-1] == LibUtilities::eQuadrilateral)
1552  {
1553  // Quad face
1554  dofs[2][faceId[i]-1] =
1555  LibUtilities::StdQuadData::getNumberOfCoefficients(P,Q) -
1556  LibUtilities::StdQuadData::getNumberOfBndCoefficients(P,Q);
1557  }
1558  else
1559  {
1560  // Tri face
1561  dofs[2][faceId[i]-1] =
1562  LibUtilities::StdTriData::getNumberOfCoefficients(P,Q) -
1563  LibUtilities::StdTriData::getNumberOfBndCoefficients(P,Q);
1564  }
1565  }
1566 
1567  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1568 
1569  // Note that nExtraDirichlet is not used in the logic below; it just
1570  // needs to be set so that the coupled solver in
1571  // IncNavierStokesSolver can work.
1572  int nExtraDirichlet;
1573  int mdswitch;
1574  m_session->LoadParameter(
1575  "MDSwitch", mdswitch, 10);
1576 
1577  int nGraphVerts =
1578  CreateGraph(locExp, bndCondExp, bndCondVec,
1579  checkIfSystemSingular, periodicVerts, periodicEdges,
1580  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1581  extraDirEdges, firstNonDirGraphVertId,
1582  nExtraDirichlet, mdswitch);
1583 
1584  /*
1585  * Set up an array which contains the offset information of the
1586  * different graph vertices.
1587  *
1588  * This basically means to identify to how many global degrees of
1589  * freedom the individual graph vertices correspond. Obviously,
1590  * the graph vertices corresponding to the mesh-vertices account
1591  * for a single global DOF. However, the graph vertices
1592  * corresponding to the element edges correspond to N-2 global DOF
1593  * where N is equal to the number of boundary modes on this edge.
1594  */
1595  Array<OneD, int> graphVertOffset(
1596  graph[0].size() + graph[1].size() + graph[2].size() + 1,0);
1597 
1598  graphVertOffset[0] = 0;
1599 
1600  for(i = 0; i < locExpVector.size(); ++i)
1601  {
1602  exp = locExpVector[i];
1603 
1604  for(j = 0; j < exp->GetNverts(); ++j)
1605  {
1606  meshVertId = exp->GetGeom()->GetVid(j);
1607  graphVertOffset[graph[0][meshVertId]+1] = 1;
1608  }
1609 
1610  for(j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1611  {
1612  if(exp->GetGeom()->GetNumFaces()) // 3D version
1613  {
1614  nEdgeInteriorCoeffs =
1615  exp->as<LocalRegions::Expansion3D>()->
1616  GetEdgeNcoeffs(j) - 2;
1617  }
1618  else
1619  {
1620  nEdgeInteriorCoeffs = exp->GetTraceNcoeffs(j) - 2;
1621  }
1622  meshEdgeId = exp->GetGeom()->GetEid(j);
1623  graphVertOffset[graph[1][meshEdgeId]+1]
1624  = dofs[1][meshEdgeId];
1625 
1626  // Need a sign vector for modal expansions if nEdgeCoeffs
1627  // >=3 (not 4 because of variable order case)
1628  if(nEdgeInteriorCoeffs &&
1629  (exp->GetBasisType(0) == LibUtilities::eModified_A))
1630  {
1631  m_signChange = true;
1632  }
1633  }
1634 
1635  for(j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1636  {
1637  meshFaceId = exp->GetGeom()->GetFid(j);
1638  graphVertOffset[graph[2][meshFaceId]+1] =
1639  dofs[2][meshFaceId];
1640  }
1641  }
1642 
1643  for(i = 1; i < graphVertOffset.size(); i++)
1644  {
1645  graphVertOffset[i] += graphVertOffset[i-1];
1646  }
1647 
1648  // Allocate the proper amount of space for the class-data
1649  m_numLocalCoeffs = numLocalCoeffs;
1650  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1651  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
1652  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
1653  m_localToLocalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
1654  m_localToLocalIntMap = Array<OneD, int>(m_numLocalCoeffs-
1655  m_numLocalBndCoeffs,-1);
1656  m_bndCondCoeffsToLocalCoeffsMap = Array<OneD, int>
1657  (m_numLocalBndCondCoeffs,-1);
1658 
1659  // If required, set up the sign-vector
1660  if(m_signChange)
1661  {
1662  m_localToGlobalSign = Array<OneD, NekDouble>
1663  (m_numLocalCoeffs,1.0);
1664  m_localToGlobalBndSign = Array<OneD, NekDouble>
1665  (m_numLocalBndCoeffs,1.0);
1666  m_bndCondCoeffsToLocalCoeffsSign = Array<OneD, NekDouble>
1667  (m_numLocalBndCondCoeffs,1.0);
1668  }
1669 
1670  m_staticCondLevel = 0;
1671  m_numPatches = locExpVector.size();
1672  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1673  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1674  for(i = 0; i < m_numPatches; ++i)
1675  {
1676  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1677  locExpVector[i]->NumBndryCoeffs();
1678  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1679  locExpVector[i]->GetNcoeffs() -
1680  locExpVector[i]->NumBndryCoeffs();
1681  }
1682 
1683  /**
1684  * STEP 6: Now, all ingredients are ready to set up the actual
1685  * local to global mapping.
1686  *
1687  * The remainder of the map consists of the element-interior
1688  * degrees of freedom. This leads to the block-diagonal submatrix
1689  * as each element-interior mode is globally orthogonal to modes
1690  * in all other elements.
1691  */
1692  cnt = 0;
1693 
1694  // Loop over all the elements in the domain
1695  int cntbdry = 0;
1696  int cntint = 0;
1697  for(i = 0; i < locExpVector.size(); ++i)
1698  {
1699  exp = locExpVector[i];
1700  cnt = locExp.GetCoeff_Offset(i);
1701 
1702  int nbdry = exp->NumBndryCoeffs();
1703  int nint = exp->GetNcoeffs() - nbdry;
1704 
1705  Array<OneD,unsigned int> bmap(nbdry);
1706  Array<OneD,unsigned int> imap(nint);
1707 
1708  exp->GetBoundaryMap(bmap);
1709  exp->GetInteriorMap(imap);
1710 
1711  for(j = 0; j < nbdry; ++j)
1712  {
1713  m_localToLocalBndMap[cntbdry++] = cnt + bmap[j];
1714  }
1715 
1716  for(j = 0; j < nint; ++j)
1717  {
1718  m_localToLocalIntMap[cntint++] = cnt + imap[j];
1719  }
1720 
1721  for(j = 0; j < exp->GetNverts(); ++j)
1722  {
1723  meshVertId = exp->GetGeom()->GetVid(j);
1724 
1725  // Set the global DOF for vertex j of element i
1726  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1727  graphVertOffset[graph[0][meshVertId]];
1728  }
1729 
1730  for(j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1731  {
1732  if(exp->GetGeom()->GetNumFaces())
1733  {
1734  nEdgeInteriorCoeffs = exp->
1735  as<LocalRegions::Expansion3D>()->
1736  GetEdgeNcoeffs(j)-2;
1737  }
1738  else
1739  {
1740  nEdgeInteriorCoeffs = exp->GetTraceNcoeffs(j)-2;
1741  }
1742  edgeOrient = exp->GetGeom()->GetEorient(j);
1743  meshEdgeId = exp->GetGeom()->GetEid(j);
1744 
1745  auto pIt = periodicEdges.find(meshEdgeId);
1746 
1747  // See if this edge is periodic. If it is, then we map all
1748  // edges to the one with lowest ID, and align all
1749  // coefficients to this edge orientation.
1750  if (pIt != periodicEdges.end())
1751  {
1752  pair<int, StdRegions::Orientation> idOrient =
1753  DeterminePeriodicEdgeOrientId(
1754  meshEdgeId, edgeOrient, pIt->second);
1755  edgeOrient = idOrient.second;
1756  }
1757 
1758  if(exp->GetGeom()->GetNumFaces())
1759  {
1760  exp->as<LocalRegions::Expansion3D>()->
1761  GetEdgeInteriorToElementMap(j,edgeInteriorMap,
1762  edgeInteriorSign,edgeOrient);
1763  }
1764  else
1765  {
1766  exp->GetTraceInteriorToElementMap(j,edgeInteriorMap,
1767  edgeInteriorSign,edgeOrient);
1768  }
1769 
1770  // Set the global DOF's for the interior modes of edge j
1771  for(k = 0; k < dofs[1][meshEdgeId]; ++k)
1772  {
1773  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1774  graphVertOffset[graph[1][meshEdgeId]]+k;
1775  }
1776  for(k = dofs[1][meshEdgeId]; k < nEdgeInteriorCoeffs; ++k)
1777  {
1778  m_localToGlobalMap[cnt+edgeInteriorMap[k]] = 0;
1779  }
1780 
1781  // Fill the sign vector if required
1782  if(m_signChange)
1783  {
1784  for(k = 0; k < dofs[1][meshEdgeId]; ++k)
1785  {
1786  m_localToGlobalSign[cnt+edgeInteriorMap[k]] =
1787  (NekDouble) edgeInteriorSign[k];
1788  }
1789  for(k = dofs[1][meshEdgeId]; k < nEdgeInteriorCoeffs; ++k)
1790  {
1791  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = 0.0;
1792  }
1793  }
1794  }
1795 
1796  for(j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1797  {
1798  faceOrient = exp->GetGeom()->GetForient(j);
1799  meshFaceId = exp->GetGeom()->GetFid(j);
1800 
1801  auto pIt = periodicFaces.find(meshFaceId);
1802 
1803  if (pIt != periodicFaces.end() &&
1804  meshFaceId == min(meshFaceId, pIt->second[0].id))
1805  {
1806  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1807  }
1808 
1809  exp->GetTraceInteriorToElementMap(j,faceInteriorMap,
1810  faceInteriorSign,
1811  faceOrient);
1812 
1813  // Set the global DOF's for the interior modes of face j
1814  exp->GetTraceNumModes(j, numModes0, numModes1, faceOrient);
1815  switch(faceType[meshFaceId])
1816  {
1817  case LibUtilities::eQuadrilateral:
1818  {
1819  int kLoc=0;
1820  k = 0;
1821  for( q = 2; q < numModes1; q++)
1822  {
1823  for( p = 2; p < numModes0; p++)
1824  {
1825  if( (p < faceModes[0][meshFaceId]) &&
1826  (q < faceModes[1][meshFaceId]))
1827  {
1828  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1829  graphVertOffset[graph[2][meshFaceId]]+k;
1830  if(m_signChange)
1831  {
1832  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1833  (NekDouble) faceInteriorSign[kLoc];
1834  }
1835  k++;
1836  }
1837  else
1838  {
1839  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] = 0;
1840  if(m_signChange)
1841  {
1842  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1843  }
1844  }
1845  kLoc++;
1846  }
1847  }
1848  }
1849  break;
1850  case LibUtilities::eTriangle:
1851  {
1852  int kLoc=0;
1853  k = 0;
1854  for( p = 2; p < numModes0; p++)
1855  {
1856  for( q = 1; q < numModes1-p; q++)
1857  {
1858  if( (p < faceModes[0][meshFaceId]) &&
1859  (p+q < faceModes[1][meshFaceId]))
1860  {
1861  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1862  graphVertOffset[graph[2][meshFaceId]]+k;
1863  if(m_signChange)
1864  {
1865  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1866  (NekDouble) faceInteriorSign[kLoc];
1867  }
1868  k++;
1869  }
1870  else
1871  {
1872  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] = 0;
1873  if(m_signChange)
1874  {
1875  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1876  }
1877  }
1878  kLoc++;
1879  }
1880  }
1881  }
1882  break;
1883  default:
1884  ASSERTL0(false,"Shape not recognised");
1885  break;
1886  }
1887  }
1888  }
1889 
1890  // Set up the mapping for the boundary conditions
1891  // Set up boundary mapping
1892  map<int, pair<int,int> > traceToElmtTraceMap;
1893  int id;
1894 
1895  for(cnt = i = 0; i < locExpVector.size(); ++i)
1896  {
1897  exp = locExpVector[i];
1898 
1899  for(j = 0; j < exp->GetNtraces(); ++j)
1900  {
1901  id = exp->GetGeom()->GetTid(j);
1902 
1903  traceToElmtTraceMap[id] = pair<int,int>(i,j);
1904  }
1905  }
1906 
1907  Array<OneD, unsigned int> maparray;
1908  Array<OneD, int> signarray;
1909  map<int,pair<int,NekDouble>> GloDirBndCoeffToLocalCoeff;
1910  set<int> CoeffOnDirTrace;
1911 
1912  cnt = 0;
1913  int offset = 0;
1914  for(i = 0; i < bndCondExp.size(); i++)
1915  {
1916  set<int> foundExtraVerts, foundExtraEdges;
1917  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1918  {
1919  bndExp = bndCondExp[i]->GetExp(j);
1920  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1921 
1922  int id = bndExp->GetGeom()->GetGlobalID();
1923 
1924  ASSERTL1(traceToElmtTraceMap.count(id) > 0,
1925  "Failed to find trace id");
1926 
1927  int eid = traceToElmtTraceMap[id].first;
1928  int tid = traceToElmtTraceMap[id].second;
1929 
1930  exp = locExpVector[eid];
1931  int dim = exp->GetShapeDimension();
1932 
1933  if(dim == 1)
1934  {
1935  m_bndCondCoeffsToLocalCoeffsMap [cnt] =
1936  locExp.GetCoeff_Offset(eid) + exp->GetVertexMap(tid);
1937  }
1938  else
1939  {
1940  if(dim == 2)
1941  {
1942  exp->GetTraceToElementMap(tid,
1943  maparray,signarray,
1944  exp->GetGeom()->
1945  GetEorient(tid),
1946  bndExp->GetBasisNumModes(0));
1947  }
1948  else if (dim == 3)
1949  {
1950  exp->GetTraceToElementMap(tid, maparray,signarray,
1951  exp->GetGeom()->
1952  GetForient(tid),
1953  bndExp->GetBasisNumModes(0),
1954  bndExp->GetBasisNumModes(1));
1955  }
1956 
1957  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1958  {
1959  m_bndCondCoeffsToLocalCoeffsMap [cnt+k] =
1960  locExp.GetCoeff_Offset(eid) + maparray[k];
1961  if(m_signChange)
1962  {
1963  m_bndCondCoeffsToLocalCoeffsSign[cnt+k] = signarray[k];
1964  }
1965  }
1966 
1967  }
1968 
1969  // we now need some information to work out how to
1970  // handle vertices and edges that are only just
1971  // touching a dirichlet boundary (and not the
1972  // whole edge/face)
1973 
1974  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1975  {
1976  int locid = m_bndCondCoeffsToLocalCoeffsMap [cnt+k];
1977  int gloid = m_localToGlobalMap[locid];
1978  NekDouble sign = 1.0;
1979 
1980  if(m_signChange)
1981  {
1982  sign = m_bndCondCoeffsToLocalCoeffsSign[cnt+k];
1983  }
1984 
1985  if (bndConditions[i]->GetBoundaryConditionType() ==
1986  SpatialDomains::eDirichlet)
1987  {
1988  CoeffOnDirTrace.insert(locid);
1989 
1990  // store the local id and sign from global id
1991  // back to local space;
1992  GloDirBndCoeffToLocalCoeff[gloid] =
1993  pair<int,NekDouble> (locid,sign);
1994  }
1995  }
1996  }
1997  offset += bndCondExp[i]->GetNcoeffs();
1998  }
1999 
2000  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
2001  m_numGlobalBndCoeffs = globalId;
2002 
2003  // Set up a mapping list of Dirichlet Local Dofs that
2004  // arise due to one vertex or edge just touching a
2005  // Dirichlet boundary and need the value from another
2006  // local coeff that has been filled by the boundary
2007  // coeffs.
2008 
2009  Array<OneD, NekDouble> gloParaDirBnd(m_numGlobalBndCoeffs,-1.0);
2010 
2011  Array<OneD, unsigned int> bndmap;
2012  cnt = 0;
2013  for(i = 0; i < locExpVector.size(); ++i)
2014  {
2015  int gloid;
2016 
2017  exp = locExpVector[i];
2018 
2019  exp->GetBoundaryMap(bndmap);
2020 
2021  for(j = 0; j < bndmap.size(); ++j)
2022  {
2023  k = cnt + bndmap[j];
2024 
2025  if(CoeffOnDirTrace.count(k) == 0)
2026  {
2027  gloid = m_localToGlobalMap[k];
2028 
2029  if(gloid < m_numGlobalDirBndCoeffs) // point on Dir BC
2030  {
2031  if(GloDirBndCoeffToLocalCoeff.count(gloid))
2032  {
2033  int locid = GloDirBndCoeffToLocalCoeff[gloid].
2034  first;
2035  NekDouble sign = 1.0;
2036 
2037  if(m_signChange)
2038  {
2039  sign = m_localToGlobalSign[locid]*
2040  m_localToGlobalSign[k];
2041  }
2042 
2043  ExtraDirDof DirDofs(k,locid,sign);
2044  // could make same `structure as extraDirDof
2045  m_copyLocalDirDofs.insert(DirDofs);
2046  }
2047  else // else could be on another parallel partition.
2048  {
2049  gloParaDirBnd[gloid] = gloid;
2050  }
2051  }
2052  }
2053  }
2054 
2055  cnt += exp->GetNcoeffs();
2056  }
2057 
2058  /*
2059  * The boundary condition mapping is generated from the same vertex
2060  * renumbering.
2061  */
2062  cnt=0;
2063  for(i = 0; i < m_numLocalCoeffs; ++i)
2064  {
2065  if(m_localToGlobalMap[i] == -1)
2066  {
2067  m_localToGlobalMap[i] = globalId++;
2068  }
2069  else
2070  {
2071  if(m_signChange)
2072  {
2073  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
2074  }
2075  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
2076  }
2077  }
2078 
2079  m_numGlobalCoeffs = globalId;
2080 
2081  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges,
2082  periodicFaces);
2083 
2084  // Now that universal map is setup reset gloParaDirBnd to
2085  // 0 if no point communicated or universal value of not
2086  // equal to -1.0
2087  for(i = 0; i < m_numGlobalBndCoeffs; ++i)
2088  {
2089  int gloid = gloParaDirBnd[i];
2090  if(gloid == -1)
2091  {
2092  gloParaDirBnd[i] = 0.0;
2093  }
2094  else
2095  {
2096  gloParaDirBnd[i] = m_globalToUniversalMap[gloid];
2097  }
2098  }
2099 
2100  // Use parallel boundary communication to set parallel
2101  // dirichlet values on all processors Needs to be after
2102  // SetupUuniversialC0ContMap
2103  Gs::Gather(gloParaDirBnd,Gs::gs_max,m_bndGsh);
2104 
2105  // copy global ids back to local values in partition to initialise gs communicator.
2106  Array<OneD, long> paraDirBnd(m_numLocalCoeffs);
2107  for(i = 0; i < numLocalCoeffs; ++i)
2108  {
2109  paraDirBnd[i] = 0.0;
2110 
2111  int id = m_localToGlobalMap[i];
2112 
2113  if(id >= m_numGlobalDirBndCoeffs)
2114  {
2115  continue;
2116  }
2117 
2118  paraDirBnd[i] = gloParaDirBnd[id];
2119 
2120  if(gloParaDirBnd[id] > 0.0)
2121  {
2122  // gather any sign changes due to edge modes
2123  if(m_signChange)
2124  {
2125  if(m_localToGlobalSign[i] < 0)
2126  {
2127  m_parallelDirBndSign.insert(i);
2128  }
2129  }
2130  }
2131  }
2132 
2133  m_dirBndGsh = Gs::Init(paraDirBnd,vComm,verbose);
2134 
2135  // Set up the local to global map for the next level when using
2136  // multi-level static condensation
2137  if ((m_solnType == eDirectMultiLevelStaticCond ||
2138  m_solnType == eIterativeMultiLevelStaticCond ||
2139  m_solnType == eXxtMultiLevelStaticCond ||
2140  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
2141  {
2142  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
2143  {
2144  Array<OneD, int> vwgts_perm(
2145  graph[0].size() + graph[1].size() + graph[2].size()
2146  - firstNonDirGraphVertId);
2147 
2148  for (i = 0; i < locExpVector.size(); ++i)
2149  {
2150  exp = locExpVector[i];
2151 
2152  for (j = 0; j < exp->GetNverts(); ++j)
2153  {
2154  meshVertId = exp->GetGeom()->GetVid(j);
2155 
2156  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
2157  {
2158  vwgts_perm[graph[0][meshVertId] -
2159  firstNonDirGraphVertId] =
2160  dofs[0][meshVertId];
2161  }
2162  }
2163 
2164  for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
2165  {
2166  meshEdgeId = exp->GetGeom()->GetEid(j);
2167 
2168  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
2169  {
2170  vwgts_perm[graph[1][meshEdgeId] -
2171  firstNonDirGraphVertId] =
2172  dofs[1][meshEdgeId];
2173  }
2174  }
2175 
2176  for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
2177  {
2178  meshFaceId = exp->GetGeom()->GetFid(j);
2179 
2180  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
2181  {
2182  vwgts_perm[graph[2][meshFaceId] -
2183  firstNonDirGraphVertId] =
2184  dofs[2][meshFaceId];
2185  }
2186  }
2187  }
2188 
2189  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
2190  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
2191  AllocateSharedPtr(this, bottomUpGraph);
2192  }
2193  }
2194 
2195  m_hash = hash_range(m_localToGlobalMap.begin(),
2196  m_localToGlobalMap.end());
2197 
2198  // Add up hash values if parallel
2199  int hash = m_hash;
2200  vComm->AllReduce(hash, LibUtilities::ReduceSum);
2201  m_hash = hash;
2202 
2203  CalculateBndSystemBandWidth();
2204  CalculateFullSystemBandWidth();
2205  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
Nektar::MultiRegions::AssemblyMapCG::CreateGraph
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)
Definition: AssemblyMapCG.cpp:80
Nektar::MultiRegions::AssemblyMapCG::m_parallelDirBndSign
std::set< int > m_parallelDirBndSign
Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change.
Definition: AssemblyMapCG.h:154
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMap
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:124
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:120
Nektar::MultiRegions::AssemblyMapCG::SetUpUniversalC0ContMap
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
Definition: AssemblyMapCG.cpp:2326
Nektar::MultiRegions::AssemblyMapCG::CalculateFullSystemBandWidth
void CalculateFullSystemBandWidth()
Calculate the bandwith of the full matrix system.
Definition: AssemblyMapCG.cpp:2724
Nektar::MultiRegions::AssemblyMapCG::m_copyLocalDirDofs
std::set< ExtraDirDof > m_copyLocalDirDofs
Set indicating degrees of freedom which are Dirichlet but whose value is stored on another processor.
Definition: AssemblyMapCG.h:151
Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:118
Nektar::MultiRegions::AssemblyMapCG::m_numLocalBndCondCoeffs
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
Definition: AssemblyMapCG.h:142
Nektar::MultiRegions::AssemblyMap::m_solnType
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:395
Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:357
Nektar::MultiRegions::AssemblyMap::m_signChange
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:371
Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToLocalCoeffsMap
Array< OneD, int > m_bndCondCoeffsToLocalCoeffsMap
Integer map of bnd cond coeffs to local coefficients.
Definition: AssemblyMap.h:382
Nektar::MultiRegions::AssemblyMap::m_numGlobalCoeffs
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:368
Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Definition: AssemblyMap.h:374
Nektar::MultiRegions::AssemblyMap::CalculateBndSystemBandWidth
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Definition: AssemblyMap.cpp:510
Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffsPerPatch
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:425
Nektar::MultiRegions::AssemblyMap::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:329
Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:338
Nektar::MultiRegions::AssemblyMap::m_nextLevelLocalToGlobalMap
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:432
Nektar::MultiRegions::AssemblyMap::m_staticCondLevel
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:421
Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:376
Nektar::MultiRegions::AssemblyMap::m_localToLocalIntMap
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
Definition: AssemblyMap.h:380
Nektar::MultiRegions::AssemblyMap::m_numGlobalDirBndCoeffs
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:344
Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:427
Nektar::MultiRegions::AssemblyMap::m_dirBndGsh
Gs::gs_data * m_dirBndGsh
gs gather communication to impose Dirhichlet BCs.
Definition: AssemblyMap.h:417
Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToLocalCoeffsSign
Array< OneD, NekDouble > m_bndCondCoeffsToLocalCoeffsSign
Integer map of sign of bnd cond coeffs to local coefficients.
Definition: AssemblyMap.h:384
Nektar::MultiRegions::AssemblyMap::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:332
Nektar::MultiRegions::AssemblyMap::m_localToLocalBndMap
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
Definition: AssemblyMap.h:378
Nektar::MultiRegions::AssemblyMap::m_numPatches
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:423
Nektar::MultiRegions::AssemblyMap::m_numGlobalBndCoeffs
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:340
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Gs::Gather
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
Gs::gs_max
@ gs_max
Definition: GsLib.hpp:53
Gs::gs_min
@ gs_min
Definition: GsLib.hpp:53
Gs::Finalise
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:226
Gs::Init
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
Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:137
Nektar::LibUtilities::StdQuadData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:147
Nektar::LibUtilities::StdTriData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:113
Nektar::LibUtilities::StdTriData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:125
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::eModifiedPyr_C
@ eModifiedPyr_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Nektar::MultiRegions::BottomUpSubStructuredGraphSharedPtr
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
Definition: SubStructuredGraph.h:59
Nektar::MultiRegions::ExtraDirDof
std::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:53
Nektar::MultiRegions::DofGraph
std::vector< std::map< int, int > > DofGraph
Definition: AssemblyMapCG.h:55
Nektar::MultiRegions::DeterminePeriodicEdgeOrientId
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...
Definition: AssemblyMapCG.cpp:2235
Nektar::MultiRegions::DeterminePeriodicFaceOrient
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:2283
Nektar::SpatialDomains::Geometry3DSharedPtr
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
Nektar::hash_range
std::size_t hash_range(Iter first, Iter last)
Definition: HashUtils.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmax
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:892
main.P
P
Definition: main.py:133
Gs::gs_data
Definition: GsLib.hpp:129

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, 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::MultiRegions::ePETScMultiLevelStaticCond, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Gs::Finalise(), Gs::Gather(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::LibUtilities::StdTriData::getNumberOfBndCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfBndCoefficients(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Gs::gs_max, Gs::gs_min, Nektar::hash_range(), Gs::Init(), Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToLocalCoeffsMap, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToLocalCoeffsSign, Nektar::MultiRegions::AssemblyMap::m_bndGsh, Nektar::MultiRegions::AssemblyMap::m_comm, m_copyLocalDirDofs, Nektar::MultiRegions::AssemblyMap::m_dirBndGsh, m_globalToUniversalMap, Nektar::MultiRegions::AssemblyMap::m_hash, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign, m_localToGlobalMap, m_localToGlobalSign, Nektar::MultiRegions::AssemblyMap::m_localToLocalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToLocalIntMap, 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, m_parallelDirBndSign, Nektar::MultiRegions::AssemblyMap::m_session, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_staticCondLevel, CellMLToNektar.cellml_metadata::p, main::P, Nektar::LibUtilities::ReduceSum, SetUpUniversalC0ContMap(), sign, and Vmath::Vmax().

◆ ~AssemblyMapCG()

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

Destructor.

Definition at line 2210 of file AssemblyMapCG.cpp.

2211  {
2212  Gs::Finalise(m_gsh);
2213  Gs::Finalise(m_bndGsh);
2214  }

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

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

2725  {
2726  int i,j;
2727  int cnt = 0;
2728  int locSize;
2729  int maxId;
2730  int minId;
2731  int bwidth = -1;
2732  for(i = 0; i < m_numPatches; ++i)
2733  {
2734  locSize = m_numLocalBndCoeffsPerPatch[i]+m_numLocalIntCoeffsPerPatch[i];
2735  maxId = -1;
2736  minId = m_numLocalCoeffs+1;
2737  for(j = 0; j < locSize; j++)
2738  {
2739  if(m_localToGlobalMap[cnt+j] >= m_numGlobalDirBndCoeffs)
2740  {
2741  if(m_localToGlobalMap[cnt+j] > maxId)
2742  {
2743  maxId = m_localToGlobalMap[cnt+j];
2744  }
2745 
2746  if(m_localToGlobalMap[cnt+j] < minId)
2747  {
2748  minId = m_localToGlobalMap[cnt+j];
2749  }
2750  }
2751  }
2752  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2753 
2754  cnt+=locSize;
2755  }
2756 
2757  m_fullSystemBandWidth = bwidth;
2758  }
Nektar::MultiRegions::AssemblyMapCG::m_fullSystemBandWidth
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
Definition: AssemblyMapCG.h:122

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

◆ 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.

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

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::StdRegions::StdExpansion::GetTraceNcoeffs(), 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().

◆ GetCopyLocalDirDofs()

std::set<ExtraDirDof>& Nektar::MultiRegions::AssemblyMapCG::GetCopyLocalDirDofs ( )
inline

Definition at line 106 of file AssemblyMapCG.h.

107  {
108  return m_copyLocalDirDofs;
109  }

References m_copyLocalDirDofs.

◆ GetParallelDirBndSign()

std::set<int>& Nektar::MultiRegions::AssemblyMapCG::GetParallelDirBndSign ( )
inline

Definition at line 111 of file AssemblyMapCG.h.

112  {
113  return m_parallelDirBndSign;
114  }

References m_parallelDirBndSign.

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

2331  {
2332  LocalRegions::ExpansionSharedPtr exp;
2333  int nVert = 0;
2334  int nEdge = 0;
2335  int nFace = 0;
2336  int maxEdgeDof = 0;
2337  int maxFaceDof = 0;
2338  int maxIntDof = 0;
2339  int dof = 0;
2340  int cnt;
2341  int i,j,k,l;
2342  int meshVertId;
2343  int meshEdgeId;
2344  int meshFaceId;
2345  int elementId;
2346  int vGlobalId;
2347  int maxBndGlobalId = 0;
2348  StdRegions::Orientation edgeOrient;
2349  StdRegions::Orientation faceOrient;
2350  Array<OneD, unsigned int> edgeInteriorMap;
2351  Array<OneD, int> edgeInteriorSign;
2352  Array<OneD, unsigned int> faceInteriorMap;
2353  Array<OneD, int> faceInteriorSign;
2354  Array<OneD, unsigned int> interiorMap;
2355 
2356  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
2357  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
2358  const bool verbose = locExp.GetSession()->DefinesCmdLineArgument("verbose");
2359 
2360  m_globalToUniversalMap = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
2361  m_globalToUniversalMapUnique = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
2362  m_globalToUniversalBndMap = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
2363  m_globalToUniversalBndMapUnique = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
2364 
2365  // Loop over all the elements in the domain to gather mesh data
2366  for(i = 0; i < locExpVector.size(); ++i)
2367  {
2368  exp = locExpVector[i];
2369 
2370  int nv = exp->GetNverts();
2371  int ne = exp->GetGeom()->GetNumEdges();
2372  int nf = exp->GetGeom()->GetNumFaces();
2373 
2374  nVert += nv;
2375  nEdge += ne;
2376  nFace += nf;
2377 
2378  // Loop over all edges (and vertices) of element i
2379  for(j = 0; j < ne; ++j)
2380  {
2381  if(nf)
2382  {
2383  dof = exp->as<LocalRegions::Expansion3D>()->
2384  GetEdgeNcoeffs(j)-2;
2385  }
2386  else
2387  {
2388  dof = exp->GetTraceNcoeffs(j)-2;
2389  }
2390 
2391  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2392  }
2393  for(j = 0; j < nf; ++j)
2394  {
2395  dof = exp->GetTraceIntNcoeffs(j);
2396  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2397  }
2398  exp->GetInteriorMap(interiorMap);
2399  dof = interiorMap.size();
2400  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2401  }
2402 
2403  // Tell other processes about how many dof we have
2404  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2405  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2406  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2407  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2408  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2409  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2410 
2411  // Assemble global to universal mapping for this process
2412  for(i = 0; i < locExpVector.size(); ++i)
2413  {
2414  exp = locExpVector[i];
2415  cnt = locExp.GetCoeff_Offset(i);
2416 
2417  int nf = exp->GetGeom()->GetNumFaces();
2418 
2419  // Loop over all vertices of element i
2420  for(j = 0; j < exp->GetNverts(); ++j)
2421  {
2422  meshVertId = exp->GetGeom()->GetVid(j);
2423  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2424 
2425  auto pIt = perVerts.find(meshVertId);
2426  if (pIt != perVerts.end())
2427  {
2428  for (k = 0; k < pIt->second.size(); ++k)
2429  {
2430  meshVertId = min(meshVertId, pIt->second[k].id);
2431  }
2432  }
2433 
2434  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2435  m_globalToUniversalBndMap[vGlobalId] =
2436  m_globalToUniversalMap[vGlobalId];
2437  maxBndGlobalId = (vGlobalId > maxBndGlobalId ?
2438  vGlobalId : maxBndGlobalId);
2439  }
2440 
2441  // Loop over all edges of element i
2442  for(j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
2443  {
2444  meshEdgeId = exp->GetGeom()->GetEid(j);
2445  auto pIt = perEdges.find(meshEdgeId);
2446  edgeOrient = exp->GetGeom()->GetEorient(j);
2447 
2448  if (pIt != perEdges.end())
2449  {
2450  pair<int, StdRegions::Orientation> idOrient =
2451  DeterminePeriodicEdgeOrientId(
2452  meshEdgeId, edgeOrient, pIt->second);
2453  meshEdgeId = idOrient.first;
2454  edgeOrient = idOrient.second;
2455  }
2456 
2457  if(nf) // 3D version
2458  {
2459  exp->as<LocalRegions::Expansion3D>()->
2460  GetEdgeInteriorToElementMap(j,edgeInteriorMap,
2461  edgeInteriorSign,edgeOrient);
2462  dof = exp->as<LocalRegions::Expansion3D>()->
2463  GetEdgeNcoeffs(j)-2;
2464  }
2465  else // 2D version
2466  {
2467  exp->GetTraceInteriorToElementMap(j,edgeInteriorMap,
2468  edgeInteriorSign,edgeOrient);
2469  dof = exp->GetTraceNcoeffs(j)-2;
2470  }
2471 
2472  // Set the global DOF's for the interior modes of edge j
2473  // for varP, ignore modes with sign == 0
2474  for(k = 0, l = 0; k < dof; ++k)
2475  {
2476  if (m_signChange)
2477  {
2478  if (m_localToGlobalSign[cnt+edgeInteriorMap[k]]==0)
2479  {
2480  continue;
2481  }
2482  }
2483  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2484  m_globalToUniversalMap[vGlobalId]
2485  = nVert + meshEdgeId * maxEdgeDof + l + 1;
2486  m_globalToUniversalBndMap[vGlobalId]=
2487  m_globalToUniversalMap[vGlobalId];
2488  maxBndGlobalId = (vGlobalId > maxBndGlobalId ?
2489  vGlobalId : maxBndGlobalId);
2490  l++;
2491  }
2492  }
2493 
2494  // Loop over all faces of element i
2495  for(j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
2496  {
2497  faceOrient = exp->GetGeom()->GetForient(j);
2498 
2499  meshFaceId = exp->GetGeom()->GetFid(j);
2500 
2501  auto pIt = perFaces.find(meshFaceId);
2502  if (pIt != perFaces.end())
2503  {
2504  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2505  {
2506  faceOrient = DeterminePeriodicFaceOrient
2507  (faceOrient,pIt->second[0].orient);
2508  }
2509  meshFaceId = min(meshFaceId, pIt->second[0].id);
2510  }
2511 
2512 
2513  exp->GetTraceInteriorToElementMap(j,faceInteriorMap,
2514  faceInteriorSign,faceOrient);
2515  dof = exp->GetTraceIntNcoeffs(j);
2516 
2517  for(k = 0, l = 0; k < dof; ++k)
2518  {
2519  if (m_signChange)
2520  {
2521  if (m_localToGlobalSign[cnt+faceInteriorMap[k]]==0)
2522  {
2523  continue;
2524  }
2525  }
2526  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2527  m_globalToUniversalMap[vGlobalId]
2528  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2529  + l + 1;
2530  m_globalToUniversalBndMap[vGlobalId]=
2531  m_globalToUniversalMap[vGlobalId];
2532 
2533  maxBndGlobalId = (vGlobalId > maxBndGlobalId ?
2534  vGlobalId : maxBndGlobalId);
2535  l++;
2536  }
2537  }
2538 
2539  // Add interior DOFs to complete universal numbering
2540  exp->GetInteriorMap(interiorMap);
2541  dof = interiorMap.size();
2542  elementId = (exp->GetGeom())->GetGlobalID();
2543  for (k = 0; k < dof; ++k)
2544  {
2545  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2546  m_globalToUniversalMap[vGlobalId]
2547  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2548  }
2549  }
2550 
2551  // Set up the GSLib universal assemble mapping
2552  // Internal DOF do not participate in any data
2553  // exchange, so we keep these set to the special GSLib id=0 so
2554  // they are ignored.
2555  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2556  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2557  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2558  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2559  {
2560  tmp[i] = m_globalToUniversalMap[i];
2561  }
2562 
2563  m_gsh = Gs::Init(tmp, vCommRow, verbose);
2564  m_bndGsh = Gs::Init(tmp2, vCommRow, verbose);
2565  Gs::Unique(tmp, vCommRow);
2566  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2567  {
2568  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2569  }
2570  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2571  {
2572  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2573  }
2574  }
Nektar::MultiRegions::AssemblyMapCG::m_globalToUniversalMapUnique
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMapCG.h:126
Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMap
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:390
Nektar::MultiRegions::AssemblyMap::m_globalToUniversalBndMapUnique
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:392
Gs::Unique
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
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436

References Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetSession(), Nektar::StdRegions::StdExpansion::GetTraceInteriorToElementMap(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), 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().

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

2888  {
2889  Array<OneD, const NekDouble> local;
2890  if(global.data() == loc.data())
2891  {
2892  local = Array<OneD, NekDouble>(m_numLocalCoeffs,loc.data());
2893  }
2894  else
2895  {
2896  local = loc; // create reference
2897  }
2898 
2899  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2900 
2901  if(m_signChange)
2902  {
2903  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2904  }
2905  else
2906  {
2907  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2908  }
2909  UniversalAssemble(global);
2910  }
Nektar::MultiRegions::AssemblyMap::UniversalAssemble
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:862
CG_Iterations.loc
loc
Definition: CG_Iterations.py:323
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:813

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

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

2915  {
2916  Assemble(loc.GetPtr(),global.GetPtr());
2917  }
Nektar::MultiRegions::AssemblyMap::Assemble
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Definition: AssemblyMap.cpp:848

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

◆ v_GetExtraDirEdges()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2981 of file AssemblyMapCG.cpp.

2982  {
2983  return m_extraDirEdges;
2984  }

References m_extraDirEdges.

◆ v_GetFullSystemBandWidth()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2941 of file AssemblyMapCG.cpp.

2942  {
2943  return m_fullSystemBandWidth;
2944  }

References m_fullSystemBandWidth.

◆ v_GetGlobalToUniversalMap() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2783 of file AssemblyMapCG.cpp.

2784  {
2785  return m_globalToUniversalMap;
2786  }

References m_globalToUniversalMap.

◆ v_GetGlobalToUniversalMap() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2766 of file AssemblyMapCG.cpp.

2767  {
2768  return m_globalToUniversalMap[i];
2769  }

References m_globalToUniversalMap.

◆ v_GetGlobalToUniversalMapUnique() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2789 of file AssemblyMapCG.cpp.

2790  {
2791  return m_globalToUniversalMapUnique;
2792  }

References m_globalToUniversalMapUnique.

◆ v_GetGlobalToUniversalMapUnique() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2771 of file AssemblyMapCG.cpp.

2772  {
2773  return m_globalToUniversalMapUnique[i];
2774  }

References m_globalToUniversalMapUnique.

◆ v_GetLocalToGlobalMap() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2777 of file AssemblyMapCG.cpp.

2778  {
2779  return m_localToGlobalMap;
2780  }

References m_localToGlobalMap.

◆ v_GetLocalToGlobalMap() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2761 of file AssemblyMapCG.cpp.

2762  {
2763  return m_localToGlobalMap[i];
2764  }

References m_localToGlobalMap.

◆ v_GetLocalToGlobalSign() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2807 of file AssemblyMapCG.cpp.

2808  {
2809  return m_localToGlobalSign;
2810  }

References m_localToGlobalSign.

◆ v_GetLocalToGlobalSign() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2794 of file AssemblyMapCG.cpp.

2796  {
2797  if(m_signChange)
2798  {
2799  return m_localToGlobalSign[i];
2800  }
2801  else
2802  {
2803  return 1.0;
2804  }
2805  }

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

◆ v_GetNumDirEdges()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2961 of file AssemblyMapCG.cpp.

2962  {
2963  return m_numDirEdges;
2964  }

References m_numDirEdges.

◆ v_GetNumDirFaces()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2966 of file AssemblyMapCG.cpp.

2967  {
2968  return m_numDirFaces;
2969  }

References m_numDirFaces.

◆ v_GetNumNonDirEdgeModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2951 of file AssemblyMapCG.cpp.

2952  {
2953  return m_numNonDirEdgeModes;
2954  }

References m_numNonDirEdgeModes.

◆ v_GetNumNonDirEdges()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2971 of file AssemblyMapCG.cpp.

2972  {
2973  return m_numNonDirEdges;
2974  }

References m_numNonDirEdges.

◆ v_GetNumNonDirFaceModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2956 of file AssemblyMapCG.cpp.

2957  {
2958  return m_numNonDirFaceModes;
2959  }

References m_numNonDirFaceModes.

◆ v_GetNumNonDirFaces()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2976 of file AssemblyMapCG.cpp.

2977  {
2978  return m_numNonDirFaces;
2979  }

References m_numNonDirFaces.

◆ v_GetNumNonDirVertexModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2946 of file AssemblyMapCG.cpp.

2947  {
2948  return m_numNonDirVertexModes;
2949  }

References m_numNonDirVertexModes.

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

2856  {
2857  Array<OneD, const NekDouble> glo;
2858  if(global.data() == loc.data())
2859  {
2860  glo = Array<OneD, NekDouble>(m_numGlobalCoeffs,global.data());
2861  }
2862  else
2863  {
2864  glo = global; // create reference
2865  }
2866 
2867 
2868  if(m_signChange)
2869  {
2870  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2871  }
2872  else
2873  {
2874  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2875  }
2876  }
Vmath::Gathr
void Gathr(int n, const T *sign, const T *x, const int *y, T *z)
Gather vector z[i] = sign[i]*x[y[i]].
Definition: Vmath.cpp:756

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

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

2881  {
2882  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2883  }
Nektar::MultiRegions::AssemblyMap::GlobalToLocal
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
Definition: AssemblyMap.cpp:834

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

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

2586  {
2587  AssemblyMapCGSharedPtr returnval;
2588 
2589  int i, j;
2590  int nverts = 0;
2591  const std::shared_ptr<LocalRegions::ExpansionVector> exp
2592  = locexp.GetExp();
2593  int nelmts = exp->size();
2594  const bool verbose = locexp.GetSession()->DefinesCmdLineArgument("verbose");
2595 
2596  // Get Default Map and turn off any searched values.
2597  returnval = MemoryManager<AssemblyMapCG>
2598  ::AllocateSharedPtr(m_session);
2599  returnval->m_solnType = solnType;
2600  returnval->m_preconType = eNull;
2601  returnval->m_maxStaticCondLevel = 0;
2602  returnval->m_signChange = false;
2603  returnval->m_comm = m_comm;
2604 
2605  // Count the number of vertices
2606  for (i = 0; i < nelmts; ++i)
2607  {
2608  nverts += (*exp)[i]->GetNverts();
2609  }
2610 
2611  returnval->m_numLocalCoeffs = nverts;
2612  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2613 
2614  // Store original global ids in this map
2615  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2616 
2617  int cnt = 0;
2618  int cnt1 = 0;
2619  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2620 
2621  // Set up local to global map;
2622  for (i = 0; i < nelmts; ++i)
2623  {
2624  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2625  {
2626  returnval->m_localToGlobalMap[cnt] =
2627  returnval->m_localToGlobalBndMap[cnt] =
2628  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2629  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2630 
2631  // Set up numLocalDirBndCoeffs
2632  if ((returnval->m_localToGlobalMap[cnt]) <
2633  m_numGlobalDirBndCoeffs)
2634  {
2635  returnval->m_numLocalDirBndCoeffs++;
2636  }
2637  cnt++;
2638  }
2639  cnt1 += (*exp)[i]->GetNcoeffs();
2640  }
2641 
2642  cnt = 0;
2643  // Reset global numbering and count number of dofs
2644  for (i = 0; i < m_numGlobalCoeffs; ++i)
2645  {
2646  if (GlobCoeffs[i] != -1)
2647  {
2648  GlobCoeffs[i] = cnt++;
2649  }
2650  }
2651 
2652  // Set up number of globalCoeffs;
2653  returnval->m_numGlobalCoeffs = cnt;
2654 
2655  // Set up number of global Dirichlet boundary coefficients
2656  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2657  {
2658  if (GlobCoeffs[i] != -1)
2659  {
2660  returnval->m_numGlobalDirBndCoeffs++;
2661  }
2662  }
2663 
2664  // Set up global to universal map
2665  if (m_globalToUniversalMap.size())
2666  {
2667  LibUtilities::CommSharedPtr vCommRow
2668  = m_session->GetComm()->GetRowComm();
2669  int nglocoeffs = returnval->m_numGlobalCoeffs;
2670  returnval->m_globalToUniversalMap
2671  = Array<OneD, int> (nglocoeffs);
2672  returnval->m_globalToUniversalMapUnique
2673  = Array<OneD, int> (nglocoeffs);
2674 
2675  // Reset local to global map and setup universal map
2676  for (i = 0; i < nverts; ++i)
2677  {
2678  cnt = returnval->m_localToGlobalMap[i];
2679  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2680 
2681  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2682  m_globalToUniversalMap[cnt];
2683  }
2684 
2685  Nektar::Array<OneD, long> tmp(nglocoeffs);
2686  Vmath::Zero(nglocoeffs, tmp, 1);
2687  for (unsigned int i = 0; i < nglocoeffs; ++i)
2688  {
2689  tmp[i] = returnval->m_globalToUniversalMap[i];
2690  }
2691  returnval->m_gsh = Gs::Init(tmp, vCommRow, verbose);
2692  Gs::Unique(tmp, vCommRow);
2693  for (unsigned int i = 0; i < nglocoeffs; ++i)
2694  {
2695  returnval->m_globalToUniversalMapUnique[i]
2696  = (tmp[i] >= 0 ? 1 : 0);
2697  }
2698  }
2699  else // not sure this option is ever needed.
2700  {
2701  for (i = 0; i < nverts; ++i)
2702  {
2703  cnt = returnval->m_localToGlobalMap[i];
2704  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2705  }
2706  }
2707 
2708  return returnval;
2709  }
Nektar::MultiRegions::AssemblyMapCGSharedPtr
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:51
Nektar::MultiRegions::eNull
@ eNull
No Solution type specified.
Definition: MultiRegions.hpp:114

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

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

2816  {
2817  Array<OneD, const NekDouble> local;
2818  if(global.data() == loc.data())
2819  {
2820  local = Array<OneD, NekDouble>(m_numLocalCoeffs,loc.data());
2821  }
2822  else
2823  {
2824  local = loc; // create reference
2825  }
2826 
2827 
2828  if(m_signChange)
2829  {
2830  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2831  }
2832  else
2833  {
2834  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2835  }
2836 
2837  // ensure all values are unique by calling a max
2838  if(useComm)
2839  {
2840  Gs::Gather(global, Gs::gs_max, m_gsh);
2841  }
2842  }
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:772

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

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

2849  {
2850  LocalToGlobal(loc.GetPtr(),global.GetPtr(),useComm);
2851  }
Nektar::MultiRegions::AssemblyMap::LocalToGlobal
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Definition: AssemblyMap.cpp:818

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

◆ v_UniversalAssemble() [1/3]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2919 of file AssemblyMapCG.cpp.

2921  {
2922  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2923  }

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

◆ v_UniversalAssemble() [2/3]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2931 of file AssemblyMapCG.cpp.

2934  {
2935  Array<OneD, NekDouble> tmp(offset);
2936  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2937  UniversalAssemble(pGlobal);
2938  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2939  }
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199

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

◆ v_UniversalAssemble() [3/3]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2925 of file AssemblyMapCG.cpp.

2927  {
2928  UniversalAssemble(pGlobal.GetPtr());
2929  }

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

Member Data Documentation

◆ m_copyLocalDirDofs

std::set<ExtraDirDof> Nektar::MultiRegions::AssemblyMapCG::m_copyLocalDirDofs
protected

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

Definition at line 151 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and GetCopyLocalDirDofs().

◆ m_extraDirEdges

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

Extra dirichlet edges in parallel.

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

Referenced by AssemblyMapCG(), 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 126 of file AssemblyMapCG.h.

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

◆ m_localToGlobalMap

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

Integer map of local coeffs to global space.

Definition at line 118 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), CalculateFullSystemBandWidth(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(), SetUpUniversalC0ContMap(), v_Assemble(), v_GetLocalToGlobalMap(), v_GlobalToLocal(), v_LinearSpaceMap(), and v_LocalToGlobal().

◆ m_localToGlobalSign

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

Integer sign of local coeffs to global space.

Definition at line 120 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(), SetUpUniversalC0ContMap(), v_Assemble(), v_GetLocalToGlobalSign(), v_GlobalToLocal(), and v_LocalToGlobal().

◆ m_maxStaticCondLevel

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

Maximum static condensation level.

Definition at line 148 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG().

◆ m_numDirEdges

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

Number of Dirichlet edges.

Definition at line 134 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 136 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 142 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 146 of file AssemblyMapCG.h.

◆ m_numNonDirEdgeModes

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

Number of non Dirichlet edge modes.

Definition at line 130 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 138 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 132 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 140 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 128 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().

◆ m_parallelDirBndSign

std::set<int> Nektar::MultiRegions::AssemblyMapCG::m_parallelDirBndSign
protected

Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change.

Definition at line 154 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and GetParallelDirBndSign().

Generated on Wed Nov 24 2021 10:33:35 for Nektar++ by doxygen 1.9.1