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 LibUtilities::CommSharedPtr &comm, 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 LibUtilities::CommSharedPtr &comm, 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 override
 
virtual int v_GetGlobalToUniversalMap (const int i) const override
 
virtual int v_GetGlobalToUniversalMapUnique (const int i) const override
 
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap () override
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap () override
 
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique () override
 
virtual NekDouble v_GetLocalToGlobalSign (const int i) const override
 
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign () const override
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const override
 
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm) const override
 
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const override
 
virtual void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const override
 
virtual void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const override
 
virtual void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const override
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const override
 
virtual void v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const override
 
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const override
 
virtual int v_GetFullSystemBandWidth () const override
 
virtual int v_GetNumNonDirVertexModes () const override
 
virtual int v_GetNumNonDirEdgeModes () const override
 
virtual int v_GetNumNonDirFaceModes () const override
 
virtual int v_GetNumDirEdges () const override
 
virtual int v_GetNumDirFaces () const override
 
virtual int v_GetNumNonDirEdges () const override
 
virtual int v_GetNumNonDirFaces () const override
 
virtual const Array< OneD, const int > & v_GetExtraDirEdges () override
 
virtual AssemblyMapSharedPtr v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType) override
 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 67 of file AssemblyMapCG.h.

Member Typedef Documentation

◆ BndCond

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

Definition at line 70 of file AssemblyMapCG.h.

◆ BndCondExp

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

Definition at line 69 of file AssemblyMapCG.h.

Constructor & Destructor Documentation

◆ AssemblyMapCG() [1/2]

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

Default constructor.

Definition at line 71 of file AssemblyMapCG.cpp.

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

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

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

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, Nektar::LibUtilities::P, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceSum, SetUpUniversalC0ContMap(), sign, and Vmath::Vmax().

◆ ~AssemblyMapCG()

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

Destructor.

Definition at line 2198 of file AssemblyMapCG.cpp.

2199 {
2200  Gs::Finalise(m_gsh);
2201  Gs::Finalise(m_bndGsh);
2202 }

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

2708 {
2709  int i, j;
2710  int cnt = 0;
2711  int locSize;
2712  int maxId;
2713  int minId;
2714  int bwidth = -1;
2715  for (i = 0; i < m_numPatches; ++i)
2716  {
2717  locSize =
2718  m_numLocalBndCoeffsPerPatch[i] + m_numLocalIntCoeffsPerPatch[i];
2719  maxId = -1;
2720  minId = m_numLocalCoeffs + 1;
2721  for (j = 0; j < locSize; j++)
2722  {
2723  if (m_localToGlobalMap[cnt + j] >= m_numGlobalDirBndCoeffs)
2724  {
2725  if (m_localToGlobalMap[cnt + j] > maxId)
2726  {
2727  maxId = m_localToGlobalMap[cnt + j];
2728  }
2729 
2730  if (m_localToGlobalMap[cnt + j] < minId)
2731  {
2732  minId = m_localToGlobalMap[cnt + j];
2733  }
2734  }
2735  }
2736  bwidth = (bwidth > (maxId - minId)) ? bwidth : (maxId - minId);
2737 
2738  cnt += locSize;
2739  }
2740 
2741  m_fullSystemBandWidth = bwidth;
2742 }
Nektar::MultiRegions::AssemblyMapCG::m_fullSystemBandWidth
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
Definition: AssemblyMapCG.h:112

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

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

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

97  {
98  return m_copyLocalDirDofs;
99  }

References m_copyLocalDirDofs.

◆ GetParallelDirBndSign()

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

Definition at line 101 of file AssemblyMapCG.h.

102  {
103  return m_parallelDirBndSign;
104  }

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

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

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2862 of file AssemblyMapCG.cpp.

2864 {
2865  Array<OneD, const NekDouble> local;
2866  if (global.data() == loc.data())
2867  {
2868  local = Array<OneD, NekDouble>(m_numLocalCoeffs, loc.data());
2869  }
2870  else
2871  {
2872  local = loc; // create reference
2873  }
2874 
2875  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2876 
2877  if (m_signChange)
2878  {
2879  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(),
2880  m_localToGlobalMap.get(), global.get());
2881  }
2882  else
2883  {
2884  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(),
2885  global.get());
2886  }
2887  UniversalAssemble(global);
2888 }
Nektar::MultiRegions::AssemblyMap::UniversalAssemble
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:806
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:862

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2890 of file AssemblyMapCG.cpp.

2892 {
2893  Assemble(loc.GetPtr(), global.GetPtr());
2894 }
Nektar::MultiRegions::AssemblyMap::Assemble
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Definition: AssemblyMap.cpp:794
Nektar::NekVector::GetPtr
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:217

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2955 of file AssemblyMapCG.cpp.

2956 {
2957  return m_extraDirEdges;
2958 }

References m_extraDirEdges.

◆ v_GetFullSystemBandWidth()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2915 of file AssemblyMapCG.cpp.

2916 {
2917  return m_fullSystemBandWidth;
2918 }

References m_fullSystemBandWidth.

◆ v_GetGlobalToUniversalMap() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2764 of file AssemblyMapCG.cpp.

2765 {
2766  return m_globalToUniversalMap;
2767 }

References m_globalToUniversalMap.

◆ v_GetGlobalToUniversalMap() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2749 of file AssemblyMapCG.cpp.

2750 {
2751  return m_globalToUniversalMap[i];
2752 }

References m_globalToUniversalMap.

◆ v_GetGlobalToUniversalMapUnique() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2769 of file AssemblyMapCG.cpp.

2771 {
2772  return m_globalToUniversalMapUnique;
2773 }

References m_globalToUniversalMapUnique.

◆ v_GetGlobalToUniversalMapUnique() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2754 of file AssemblyMapCG.cpp.

2755 {
2756  return m_globalToUniversalMapUnique[i];
2757 }

References m_globalToUniversalMapUnique.

◆ v_GetLocalToGlobalMap() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2759 of file AssemblyMapCG.cpp.

2760 {
2761  return m_localToGlobalMap;
2762 }

References m_localToGlobalMap.

◆ v_GetLocalToGlobalMap() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2744 of file AssemblyMapCG.cpp.

2745 {
2746  return m_localToGlobalMap[i];
2747 }

References m_localToGlobalMap.

◆ v_GetLocalToGlobalSign() [1/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2787 of file AssemblyMapCG.cpp.

2788 {
2789  return m_localToGlobalSign;
2790 }

References m_localToGlobalSign.

◆ v_GetLocalToGlobalSign() [2/2]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2775 of file AssemblyMapCG.cpp.

2776 {
2777  if (m_signChange)
2778  {
2779  return m_localToGlobalSign[i];
2780  }
2781  else
2782  {
2783  return 1.0;
2784  }
2785 }

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

◆ v_GetNumDirEdges()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2935 of file AssemblyMapCG.cpp.

2936 {
2937  return m_numDirEdges;
2938 }

References m_numDirEdges.

◆ v_GetNumDirFaces()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2940 of file AssemblyMapCG.cpp.

2941 {
2942  return m_numDirFaces;
2943 }

References m_numDirFaces.

◆ v_GetNumNonDirEdgeModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2925 of file AssemblyMapCG.cpp.

2926 {
2927  return m_numNonDirEdgeModes;
2928 }

References m_numNonDirEdgeModes.

◆ v_GetNumNonDirEdges()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2945 of file AssemblyMapCG.cpp.

2946 {
2947  return m_numNonDirEdges;
2948 }

References m_numNonDirEdges.

◆ v_GetNumNonDirFaceModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2930 of file AssemblyMapCG.cpp.

2931 {
2932  return m_numNonDirFaceModes;
2933 }

References m_numNonDirFaceModes.

◆ v_GetNumNonDirFaces()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2950 of file AssemblyMapCG.cpp.

2951 {
2952  return m_numNonDirFaces;
2953 }

References m_numNonDirFaces.

◆ v_GetNumNonDirVertexModes()

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2920 of file AssemblyMapCG.cpp.

2921 {
2922  return m_numNonDirVertexModes;
2923 }

References m_numNonDirVertexModes.

◆ v_GlobalToLocal() [1/2]

void Nektar::MultiRegions::AssemblyMapCG::v_GlobalToLocal ( const Array< OneD, const NekDouble > &  global,
Array< OneD, NekDouble > &  loc 
) const
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2831 of file AssemblyMapCG.cpp.

2833 {
2834  Array<OneD, const NekDouble> glo;
2835  if (global.data() == loc.data())
2836  {
2837  glo = Array<OneD, NekDouble>(m_numGlobalCoeffs, global.data());
2838  }
2839  else
2840  {
2841  glo = global; // create reference
2842  }
2843 
2844  if (m_signChange)
2845  {
2846  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(),
2847  m_localToGlobalMap.get(), loc.get());
2848  }
2849  else
2850  {
2851  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(),
2852  loc.get());
2853  }
2854 }
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:805

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2856 of file AssemblyMapCG.cpp.

2858 {
2859  GlobalToLocal(global.GetPtr(), loc.GetPtr());
2860 }
Nektar::MultiRegions::AssemblyMap::GlobalToLocal
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
Definition: AssemblyMap.cpp:782

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 
)
overrideprotectedvirtual

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

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

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2792 of file AssemblyMapCG.cpp.

2795 {
2796  Array<OneD, const NekDouble> local;
2797  if (global.data() == loc.data())
2798  {
2799  local = Array<OneD, NekDouble>(m_numLocalCoeffs, loc.data());
2800  }
2801  else
2802  {
2803  local = loc; // create reference
2804  }
2805 
2806  if (m_signChange)
2807  {
2808  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(),
2809  m_localToGlobalMap.get(), global.get());
2810  }
2811  else
2812  {
2813  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(),
2814  global.get());
2815  }
2816 
2817  // ensure all values are unique by calling a max
2818  if (useComm)
2819  {
2820  Gs::Gather(global, Gs::gs_max, m_gsh);
2821  }
2822 }
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:822

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2824 of file AssemblyMapCG.cpp.

2827 {
2828  LocalToGlobal(loc.GetPtr(), global.GetPtr(), useComm);
2829 }
Nektar::MultiRegions::AssemblyMap::LocalToGlobal
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Definition: AssemblyMap.cpp:768

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2896 of file AssemblyMapCG.cpp.

2897 {
2898  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2899 }

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
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2906 of file AssemblyMapCG.cpp.

2908 {
2909  Array<OneD, NekDouble> tmp(offset);
2910  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2911  UniversalAssemble(pGlobal);
2912  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2913 }
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

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

◆ v_UniversalAssemble() [3/3]

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

Reimplemented from Nektar::MultiRegions::AssemblyMap.

Definition at line 2901 of file AssemblyMapCG.cpp.

2902 {
2903  UniversalAssemble(pGlobal.GetPtr());
2904 }

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

Referenced by AssemblyMapCG().

◆ m_numDirEdges

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

Number of Dirichlet edges.

Definition at line 124 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirEdges().

◆ m_numDirFaces

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

Number of Dirichlet faces.

Definition at line 126 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumDirFaces().

◆ m_numLocalBndCondCoeffs

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

Number of local boundary condition coefficients.

Definition at line 132 of file AssemblyMapCG.h.

Referenced by AssemblyMapCG(), and CreateGraph().

◆ m_numLocDirBndCondDofs

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

Number of local boundary condition degrees of freedom.

Definition at line 136 of file AssemblyMapCG.h.

◆ m_numNonDirEdgeModes

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

Number of non Dirichlet edge modes.

Definition at line 120 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().

◆ m_numNonDirEdges

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

Number of Dirichlet edges.

Definition at line 128 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirEdges().

◆ m_numNonDirFaceModes

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

Number of non Dirichlet face modes.

Definition at line 122 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().

◆ m_numNonDirFaces

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

Number of Dirichlet faces.

Definition at line 130 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirFaces().

◆ m_numNonDirVertexModes

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

Number of non Dirichlet vertex modes.

Definition at line 118 of file AssemblyMapCG.h.

Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().

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

Referenced by AssemblyMapCG(), and GetParallelDirBndSign().

Generated on Sat Mar 25 2023 12:35:59 for Nektar++ by doxygen 1.9.1