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

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

#include <AssemblyMapCG2D.h>

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

Public Member Functions

 AssemblyMapCG2D (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Default constructor.
 AssemblyMapCG2D (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const Array< OneD, const ExpListSharedPtr > &bndCondExp, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions, const PeriodicMap &periodicVertsId, const PeriodicMap &periodicEdgesId, const bool checkIfSystemSingular, const std::string variable="DefaultVar")
 Constructor for the 2D expansion mappings with boundary conditions.
 AssemblyMapCG2D (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const std::string variable="DefaultVar")
 General constructor for expansions of all dimensions without boundary conditions.
virtual ~AssemblyMapCG2D ()
 Destructor.
int SetUp2DGraphC0ContMap (const ExpList &locExp, const Array< OneD, const MultiRegions::ExpListSharedPtr > &bndCondExp, const Array< OneD, Array< OneD, const SpatialDomains::BoundaryConditionShPtr > > &bndConditions, const PeriodicMap &periodicVertsId, const PeriodicMap &periodicEdgesId, Array< OneD, map< int, int > > &Dofs, Array< OneD, map< int, int > > &ReorderedGraphVertId, int &firstNonDirGraphVertID, int &nExtraDirichlet, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, set< int > &extraDirVerts, const bool checkIfSystemSingular=false, int mdswitch=1, bool doInteriorMap=false)
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Default constructor.
 AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp)
 General constructor for expansions of all dimensions without boundary conditions.
virtual ~AssemblyMapCG ()
 Destructor.
map< int, vector< ExtraDirDof > > & GetExtraDirDofs ()
- Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
 AssemblyMap ()
 Default constructor.
 AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
 Constructor with a communicator.
 AssemblyMap (AssemblyMap *oldLevelMap, const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
 Constructor for next level in multi-level static condensation.
virtual ~AssemblyMap ()
 Destructor.
LibUtilities::CommSharedPtr GetComm ()
 Retrieves the communicator.
size_t GetHash () const
 Retrieves the hash of this map.
int GetLocalToGlobalMap (const int i) const
int GetGlobalToUniversalMap (const int i) const
int GetGlobalToUniversalMapUnique (const int i) const
const Array< OneD, const int > & GetLocalToGlobalMap ()
const Array< OneD, const int > & GetGlobalToUniversalMap ()
const Array< OneD, const int > & GetGlobalToUniversalMapUnique ()
NekDouble GetLocalToGlobalSign (const int i) const
const Array< OneD, NekDouble > & GetLocalToGlobalSign () const
void LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
void LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
void GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
void GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
void Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
void Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
void UniversalAssemble (NekVector< NekDouble > &pGlobal) const
void UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
int GetLocalToGlobalBndMap (const int i) const
 Retrieve the global index of a given local boundary mode.
const Array< OneD, const int > & GetLocalToGlobalBndMap ()
 Retrieve the global indices of the local boundary modes.
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.
NekDouble GetLocalToGlobalBndSign (const int i) const
 Retrieve the sign change of a given local boundary mode.
Array< OneD, const NekDouble > GetLocalToGlobalBndSign () const
 Retrieve the sign change for all local boundary modes.
int GetBndCondCoeffsToGlobalCoeffsMap (const int i)
 Retrieves the global index corresponding to a boundary expansion mode.
const Array< OneD, const int > & GetBndCondCoeffsToGlobalCoeffsMap ()
 Retrieves the global indices corresponding to the boundary expansion modes.
NekDouble GetBndCondCoeffsToGlobalCoeffsSign (const int i)
 Returns the modal sign associated with a given boundary expansion mode.
int GetBndCondTraceToGlobalTraceMap (const int i)
 Returns the global index of the boundary trace giving the index on the boundary expansion.
const Array< OneD, const int > & GetBndCondTraceToGlobalTraceMap ()
int GetNumGlobalDirBndCoeffs () const
 Returns the number of global Dirichlet boundary coefficients.
int GetNumLocalDirBndCoeffs () const
 Returns the number of local Dirichlet boundary coefficients.
int GetNumGlobalBndCoeffs () const
 Returns the total number of global boundary coefficients.
int GetNumLocalBndCoeffs () const
 Returns the total number of local boundary coefficients.
int GetNumLocalCoeffs () const
 Returns the total number of local coefficients.
int GetNumGlobalCoeffs () const
 Returns the total number of global coefficients.
bool GetSingularSystem () const
 Retrieves if the system is singular (true) or not (false)
void GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const
void GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
void GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc, int offset) const
void GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
void LocalBndToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
void LocalBndToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const
void LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
void AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
void AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
void AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const
void AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
void UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal) const
void UniversalAssembleBnd (NekVector< NekDouble > &pGlobal) const
void UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal, int offset) const
int GetFullSystemBandWidth () const
int GetNumNonDirVertexModes () const
int GetNumNonDirEdgeModes () const
int GetNumNonDirFaceModes () const
int GetNumDirEdges () const
int GetNumDirFaces () const
int GetNumNonDirEdges () const
int GetNumNonDirFaces () const
void PrintStats (std::ostream &out, std::string variable) const
const Array< OneD, const int > & GetExtraDirEdges ()
boost::shared_ptr< AssemblyMap > LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
int GetBndSystemBandWidth () const
 Returns the bandwidth of the boundary system.
int GetStaticCondLevel () const
 Returns the level of static condensation for this map.
int GetNumPatches () const
 Returns the number of patches in this static condensation level.
const Array< OneD, const
unsigned int > & 		GetNumLocalBndCoeffsPerPatch ()
 Returns the number of local boundary coefficients in each patch.
const Array< OneD, const
unsigned int > & 		GetNumLocalIntCoeffsPerPatch ()
 Returns the number of local interior coefficients in each patch.
const AssemblyMapSharedPtr GetNextLevelLocalToGlobalMap () const
 Returns the local to global mapping for the next level in the multi-level static condensation.
void SetNextLevelLocalToGlobalMap (AssemblyMapSharedPtr pNextLevelLocalToGlobalMap)
const PatchMapSharedPtr & GetPatchMapFromPrevLevel (void) const
 Returns the patch map from the previous level of the multi-level static condensation.
bool AtLastLevel () const
 Returns true if this is the last level in the multi-level static condensation.
GlobalSysSolnType GetGlobalSysSolnType () const
 Returns the method of solving global systems.
PreconditionerType GetPreconType () const
NekDouble GetIterativeTolerance () const
int GetSuccessiveRHS () const
int GetLowestStaticCondLevel () const

Private Member Functions

void SetUp2DExpansionC0ContMap (const int numLocalCoeffs, const ExpList &locExp, const Array< OneD, const MultiRegions::ExpListSharedPtr > &bndCondExp=NullExpListSharedPtrArray, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions=SpatialDomains::NullBoundaryConditionShPtrArray, const PeriodicMap &periodicVertsId=NullPeriodicMap, const PeriodicMap &periodicEdgesId=NullPeriodicMap, const bool checkIfSystemSingular=false)
 Construct mappings for a two-dimensional scalar expansion.

Additional Inherited Members

- Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
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.
virtual int v_GetLocalToGlobalMap (const int i) const
virtual int v_GetGlobalToUniversalMap (const int i) const
virtual int v_GetGlobalToUniversalMapUnique (const int i) const
virtual const Array< OneD,
const int > & 		v_GetLocalToGlobalMap ()
virtual const Array< OneD,
const int > & 		v_GetGlobalToUniversalMap ()
virtual const Array< OneD,
const int > & 		v_GetGlobalToUniversalMapUnique ()
virtual NekDouble v_GetLocalToGlobalSign (const int i) const
virtual const Array< OneD,
NekDouble > & 		v_GetLocalToGlobalSign () const
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
virtual void v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
virtual void v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const
virtual void v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
virtual void v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const
virtual void v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const
virtual void v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const
virtual int v_GetFullSystemBandWidth () const
virtual int v_GetNumNonDirVertexModes () const
virtual int v_GetNumNonDirEdgeModes () const
virtual int v_GetNumNonDirFaceModes () const
virtual int v_GetNumDirEdges () const
virtual int v_GetNumDirFaces () const
virtual int v_GetNumNonDirEdges () const
virtual int v_GetNumNonDirFaces () const
virtual const Array< OneD,
const int > & 		v_GetExtraDirEdges ()
virtual AssemblyMapSharedPtr v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
 Construct an AssemblyMapCG object which corresponds to the linear space of the current object.
- Protected Attributes inherited from Nektar::MultiRegions::AssemblyMapCG
Array< OneD, int > m_localToGlobalMap
 Integer map of local coeffs to global space.
Array< OneD, NekDouble > m_localToGlobalSign
 Integer sign of local coeffs to global space.
int m_fullSystemBandWidth
 Bandwith of the full matrix system (no static condensation).
Array< OneD, int > m_globalToUniversalMap
 Integer map of process coeffs to universal space.
Array< OneD, int > m_globalToUniversalMapUnique
 Integer map of unique process coeffs to universal space (signed)
int m_numNonDirVertexModes
 Number of non Dirichlet vertex modes.
int m_numNonDirEdgeModes
 Number of non Dirichlet edge modes.
int m_numNonDirFaceModes
 Number of non Dirichlet face modes.
int m_numDirEdges
 Number of Dirichlet edges.
int m_numDirFaces
 Number of Dirichlet faces.
int m_numNonDirEdges
 Number of Dirichlet edges.
int m_numNonDirFaces
 Number of Dirichlet faces.
Array< OneD, int > m_extraDirEdges
 Extra dirichlet edges in parallel.
int m_maxStaticCondLevel
 Maximum static condensation level.
map< int, vector< ExtraDirDof > > m_extraDirDofs
 Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor.

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 55 of file AssemblyMapCG2D.h.

Constructor & Destructor Documentation

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

Default constructor.

Definition at line 70 of file AssemblyMapCG2D.cpp.

:
AssemblyMapCG(pSession,variable)
{
}
Nektar::MultiRegions::AssemblyMapCG2D::AssemblyMapCG2D ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const int  numLocalCoeffs,
const ExpList &  locExp,
const Array< OneD, const ExpListSharedPtr > &  bndCondExp,
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &  bndConditions,
const PeriodicMap &  periodicVertsId,
const PeriodicMap &  periodicEdgesId,
const bool  checkIfSystemSingular,
const std::string  variable = "DefaultVar" 
)

Constructor for the 2D expansion mappings with boundary conditions.

Definition at line 83 of file AssemblyMapCG2D.cpp.

References Nektar::MultiRegions::AssemblyMap::CalculateBndSystemBandWidth(), Nektar::MultiRegions::AssemblyMapCG::CalculateFullSystemBandWidth(), and SetUp2DExpansionC0ContMap().

:
AssemblyMapCG(pSession,variable)
{
SetUp2DExpansionC0ContMap(numLocalCoeffs,
locExp,
bndCondExp,
bndConditions,
periodicVertsId,
periodicEdgesId,
checkIfSystemSingular);
CalculateBndSystemBandWidth();
CalculateFullSystemBandWidth();
}
Nektar::MultiRegions::AssemblyMapCG2D::AssemblyMapCG2D ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const int  numLocalCoeffs,
const ExpList &  locExp,
const std::string  variable = "DefaultVar" 
)

General constructor for expansions of all dimensions without boundary conditions.

Definition at line 112 of file AssemblyMapCG2D.cpp.

References Nektar::MultiRegions::AssemblyMap::CalculateBndSystemBandWidth(), Nektar::MultiRegions::AssemblyMapCG::CalculateFullSystemBandWidth(), and SetUp2DExpansionC0ContMap().

:
AssemblyMapCG(pSession,variable)
{
SetUp2DExpansionC0ContMap(numLocalCoeffs, locExp);
CalculateBndSystemBandWidth();
CalculateFullSystemBandWidth();
}
Nektar::MultiRegions::AssemblyMapCG2D::~AssemblyMapCG2D ( )
virtual

Destructor.

Definition at line 128 of file AssemblyMapCG2D.cpp.

{
}

Member Function Documentation

void Nektar::MultiRegions::AssemblyMapCG2D::SetUp2DExpansionC0ContMap ( const int  numLocalCoeffs,
const ExpList &  locExp,
const Array< OneD, const MultiRegions::ExpListSharedPtr > &  bndCondExp = NullExpListSharedPtrArray,
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &  bndConditions = SpatialDomains::NullBoundaryConditionShPtrArray,
const PeriodicMap &  periodicVertsId = NullPeriodicMap,
const PeriodicMap &  periodicEdgesId = NullPeriodicMap,
const bool  checkIfSystemSingular = false 
)
private

Construct mappings for a two-dimensional scalar expansion.

Construction of the local->global map is achieved in several stages. A mesh vertex and mesh edge renumbering is constructed in #vertReorderedGraphVertId and #edgeReorderedGraphVertId through a call ot SetUp2DGraphC0ContMap.

The local numbering is then deduced in the following steps:

STEP 1: Wrap boundary conditions in an array (since routine is set up for multiple fields) and call the graph re-ordering subroutine to obtain the reordered values

STEP 2: Count out the number of Dirichlet vertices and edges first

STEP 3: Set up an array which contains the offset information of the different graph vertices.

This basically means to identify to how many global degrees of freedom the individual graph vertices correspond. Obviously, the graph vertices corresponding to the mesh-vertices account for a single global DOF. However, the graph vertices corresponding to the element edges correspond to N-2 global DOF where N is equal to the number of boundary modes on this edge.

STEP 4: 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.

STEP 5: The boundary condition mapping is generated from the same vertex renumbering.

Definition at line 142 of file AssemblyMapCG2D.cpp.

References Nektar::StdRegions::eBackwards, Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::StdRegions::eForwards, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), Nektar::StdRegions::StdExpansion::GetVertexMap(), Nektar::iterator, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsMap, Nektar::MultiRegions::AssemblyMap::m_bndCondCoeffsToGlobalCoeffsSign, Nektar::MultiRegions::AssemblyMap::m_comm, Nektar::MultiRegions::AssemblyMapCG::m_extraDirDofs, Nektar::MultiRegions::AssemblyMap::m_hash, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndMap, Nektar::MultiRegions::AssemblyMap::m_localToGlobalBndSign, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalMap, Nektar::MultiRegions::AssemblyMapCG::m_localToGlobalSign, Nektar::MultiRegions::AssemblyMapCG::m_maxStaticCondLevel, 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, Nektar::MultiRegions::AssemblyMap::m_numLocalCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalDirBndCoeffs, Nektar::MultiRegions::AssemblyMap::m_numLocalIntCoeffsPerPatch, Nektar::MultiRegions::AssemblyMap::m_numPatches, Nektar::MultiRegions::AssemblyMap::m_signChange, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::AssemblyMap::m_staticCondLevel, Nektar::MultiRegions::AssemblyMap::m_systemSingular, Nektar::NullNekDouble1DArray, Nektar::LibUtilities::ReduceSum, SetUp2DGraphC0ContMap(), Nektar::MultiRegions::AssemblyMapCG::SetUpUniversalC0ContMap(), Nektar::MultiRegions::AssemblyMap::UniversalAssembleBnd(), and Vmath::Vmax().

Referenced by AssemblyMapCG2D().

{
int i,j,k;
int cnt = 0,offset=0;
int meshVertId;
int meshEdgeId;
int bndEdgeCnt;
int globalId, nGraphVerts;
int nEdgeCoeffs;
int nEdgeInteriorCoeffs;
int firstNonDirGraphVertId;
int nLocBndCondDofs = 0;
int nLocDirBndCondDofs = 0;
int nExtraDirichlet = 0;
LocalRegions::Expansion2DSharedPtr locExpansion;
LocalRegions::SegExpSharedPtr bndSegExp;
LibUtilities::BasisType bType;
StdRegions::Orientation edgeOrient;
Array<OneD, unsigned int> edgeInteriorMap;
Array<OneD, int> edgeInteriorSign;
const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
m_signChange = false;
m_systemSingular = false;
Array<OneD, map<int,int> > ReorderedGraphVertId(2);
Array<OneD, map<int,int> > Dofs(2);
BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
set<int> extraDirVerts;
LocalRegions::Expansion2DSharedPtr exp2d;
for(i = 0; i < locExpVector.size(); ++i)
{
exp2d = locExpVector[i]->as<LocalRegions::Expansion2D>();
for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
{
// Vert Dofs
Dofs[0][(exp2d->GetGeom2D())->GetVid(j)] = 1;
// Edge Dofs
Dofs[1][(exp2d->GetGeom2D())->GetEid(j)] =
exp2d->GetEdgeNcoeffs(j)-2;
}
}
/**
* STEP 1: Wrap boundary conditions in an array (since
* routine is set up for multiple fields) and call the
* graph re-ordering subroutine to obtain the reordered
* values
*/
Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> > bndConditionsVec(1,bndConditions);
nGraphVerts = SetUp2DGraphC0ContMap(locExp,
bndCondExp,bndConditionsVec,
periodicVertsId,
periodicEdgesId,
Dofs,
ReorderedGraphVertId,
firstNonDirGraphVertId,
nExtraDirichlet,
bottomUpGraph,
extraDirVerts,
checkIfSystemSingular);
/**
* STEP 2: Count out the number of Dirichlet vertices and
* edges first
*/
for(i = 0; i < bndCondExp.num_elements(); i++)
{
for(j = 0; j < bndCondExp[i]->GetExpSize(); j++)
{
bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
if(bndConditions[i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
{
nLocDirBndCondDofs += bndSegExp->GetNcoeffs();
}
if( bndCondExp[i]->GetExp(j)==bndCondExp[bndCondExp.num_elements()-1]->GetExp(bndCondExp[bndCondExp.num_elements()-1]->GetExpSize()-1)
&& i==(bndCondExp.num_elements()-1)
&& nLocDirBndCondDofs ==0
&& checkIfSystemSingular==true)
{
nLocDirBndCondDofs =1;
m_systemSingular=true;
}
nLocBndCondDofs += bndSegExp->GetNcoeffs();
}
}
m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
/**
* STEP 3: Set up an array which contains the offset information of
* the different graph vertices.
*
* This basically means to identify to how many global degrees of
* freedom the individual graph vertices correspond. Obviously,
* the graph vertices corresponding to the mesh-vertices account
* for a single global DOF. However, the graph vertices
* corresponding to the element edges correspond to N-2 global DOF
* where N is equal to the number of boundary modes on this edge.
*/
Array<OneD, int> graphVertOffset(ReorderedGraphVertId[0].size()+
ReorderedGraphVertId[1].size()+1);
graphVertOffset[0] = 0;
m_numLocalBndCoeffs = 0;
for(i = 0; i < locExpVector.size(); ++i)
{
locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]->as<LocalRegions::Expansion2D>();
for(j = 0; j < locExpansion->GetNedges(); ++j)
{
nEdgeCoeffs = locExpansion->GetEdgeNcoeffs(j);
meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
graphVertOffset[ReorderedGraphVertId[0][meshVertId]+1] = 1;
graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]+1] = nEdgeCoeffs-2;
bType = locExpansion->GetEdgeBasisType(j);
// need a sign vector for modal expansions if nEdgeCoeffs >=4
if( (nEdgeCoeffs >= 4)&&
( (bType == LibUtilities::eModified_A)||
(bType == LibUtilities::eModified_B) ) )
{
m_signChange = true;
}
}
m_numLocalBndCoeffs += locExpVector[i]->NumBndryCoeffs();
}
for(i = 1; i < graphVertOffset.num_elements(); i++)
{
graphVertOffset[i] += graphVertOffset[i-1];
}
// Allocate the proper amount of space for the class-data and fill
// information that is already known
m_numLocalCoeffs = numLocalCoeffs;
m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(nLocBndCondDofs,-1);
// If required, set up the sign-vector
if(m_signChange)
{
m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD, NekDouble>(nLocBndCondDofs,1.0);
}
else
{
m_localToGlobalSign = NullNekDouble1DArray;
m_localToGlobalBndSign = NullNekDouble1DArray;
m_bndCondCoeffsToGlobalCoeffsSign = NullNekDouble1DArray;
}
// Set up information for multi-level static condensation.
m_staticCondLevel = 0;
m_numPatches = locExpVector.size();
m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
for(i = 0; i < m_numPatches; ++i)
{
int elmtid = locExp.GetOffset_Elmt_Id(i);
locExpansion = locExpVector[elmtid]->as<LocalRegions::Expansion2D>();
m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
locExpVector[elmtid]->NumBndryCoeffs();
m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
locExpVector[elmtid]->GetNcoeffs() -
locExpVector[elmtid]->NumBndryCoeffs();
}
/**
* STEP 4: 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.
*/
cnt = 0;
// Loop over all the elements in the domain
for(i = 0; i < locExpVector.size(); ++i)
{
locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
cnt = locExp.GetCoeff_Offset(i);
// Loop over all edges (and vertices) of element i
for(j = 0; j < locExpansion->GetNedges(); ++j)
{
PeriodicMap::const_iterator it;
nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
edgeOrient = (locExpansion->GetGeom2D())->GetEorient(j);
meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
/*
* Where the edge orientations of periodic edges matches,
* we reverse the vertices of each edge in
* DisContField2D::GetPeriodicEdges. We must therefore
* reverse the orientation of precisely one of the two
* edges so that the sign array is correctly populated.
*/
it = periodicEdgesId.find(meshEdgeId);
if (it != periodicEdgesId.end() &&
it->second[0].orient == StdRegions::eBackwards &&
meshEdgeId == min(meshEdgeId, it->second[0].id))
{
edgeOrient = edgeOrient == StdRegions::eForwards ?
StdRegions::eBackwards :
StdRegions::eForwards;
}
locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
// Set the global DOF for vertex j of element i
m_localToGlobalMap[cnt+locExpansion->GetVertexMap(j)] =
graphVertOffset[ReorderedGraphVertId[0][meshVertId]];
// Set the global DOF's for the interior modes of edge j
for(k = 0; k < nEdgeInteriorCoeffs; ++k)
{
m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]]+k;
}
// Fill the sign vector if required
if(m_signChange)
{
for(k = 0; k < nEdgeInteriorCoeffs; ++k)
{
m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
}
}
}
cnt += locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs();
}
// Set up the mapping for the boundary conditions
offset = cnt = 0;
for(i = 0; i < bndCondExp.num_elements(); i++)
{
set<int> foundExtraVerts;
for(j = 0; j < bndCondExp[i]->GetExpSize(); j++)
{
bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
for(k = 0; k < 2; k++)
{
meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]];
set<int>::iterator iter = extraDirVerts.find(meshVertId);
if (iter != extraDirVerts.end() &&
foundExtraVerts.count(meshVertId) == 0)
{
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
bndSegExp->GetVertexMap(k);
int gid = graphVertOffset[
ReorderedGraphVertId[0][meshVertId]];
ExtraDirDof t(loc, gid, 1.0);
m_extraDirDofs[i].push_back(t);
foundExtraVerts.insert(meshVertId);
}
}
meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
bndEdgeCnt = 0;
for(k = 0; k < bndSegExp->GetNcoeffs(); k++)
{
if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
{
m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]]+bndEdgeCnt;
bndEdgeCnt++;
}
}
}
offset += bndCondExp[i]->GetNcoeffs();
}
globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
m_numGlobalBndCoeffs = globalId;
/**
* STEP 5: The boundary condition mapping is generated from the
* same vertex renumbering.
*/
cnt=0;
for(i = 0; i < m_numLocalCoeffs; ++i)
{
if(m_localToGlobalMap[i] == -1)
{
m_localToGlobalMap[i] = globalId++;
}
else
{
if(m_signChange)
{
m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
}
m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
}
}
m_numGlobalCoeffs = globalId;
SetUpUniversalC0ContMap(locExp, periodicVertsId, periodicEdgesId);
// Since we can now have multiple entries to m_extraDirDofs due to
// periodic boundary conditions we make a call to work out the
// multiplicity of all entries and invert
Array<OneD, NekDouble> valence(m_numGlobalBndCoeffs,0.0);
// Fill in Dirichlet coefficients that are to be sent to other
// processors with a value of 1
map<int, vector<ExtraDirDof> >::iterator Tit;
// Generate valence for extraDirDofs
for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
{
for (i = 0; i < Tit->second.size(); ++i)
{
valence[Tit->second[i].get<1>()] = 1.0;
}
}
UniversalAssembleBnd(valence);
// Set third argument of tuple to inverse of valence.
for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
{
for (i = 0; i < Tit->second.size(); ++i)
{
boost::get<2>(Tit->second.at(i)) = 1.0/valence[Tit->second.at(i).get<1>()];
}
}
// Set up the local to global map for the next level when using
// multi-level static condensation
if ((m_solnType == eDirectMultiLevelStaticCond ||
m_solnType == eIterativeMultiLevelStaticCond ||
m_solnType == eXxtMultiLevelStaticCond) && nGraphVerts)
{
if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1) &&
m_staticCondLevel < m_maxStaticCondLevel)
{
Array<OneD, int> vwgts_perm(
Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
for(i = 0; i < locExpVector.size(); ++i)
{
int eid = locExp.GetOffset_Elmt_Id(i);
locExpansion = locExpVector[eid]->as<LocalRegions::Expansion2D>();
for(j = 0; j < locExpansion->GetNverts(); ++j)
{
meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
if(ReorderedGraphVertId[0][meshVertId] >=
firstNonDirGraphVertId)
{
vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
firstNonDirGraphVertId] =
Dofs[0][meshVertId];
}
if(ReorderedGraphVertId[1][meshEdgeId] >=
firstNonDirGraphVertId)
{
vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
firstNonDirGraphVertId] =
Dofs[1][meshEdgeId];
}
}
}
bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
AllocateSharedPtr(this,bottomUpGraph);
}
}
m_hash = boost::hash_range(
m_localToGlobalMap.begin(), m_localToGlobalMap.end());
// Add up hash values if parallel
int hash = m_hash;
m_comm->GetRowComm()->AllReduce(hash,
LibUtilities::ReduceSum);
m_hash = hash;
}
int Nektar::MultiRegions::AssemblyMapCG2D::SetUp2DGraphC0ContMap ( const ExpList &  locExp,
const Array< OneD, const MultiRegions::ExpListSharedPtr > &  bndCondExp,
const Array< OneD, Array< OneD, const SpatialDomains::BoundaryConditionShPtr > > &  bndConditions,
const PeriodicMap &  periodicVertsId,
const PeriodicMap &  periodicEdgesId,
Array< OneD, map< int, int > > &  Dofs,
Array< OneD, map< int, int > > &  ReorderedGraphVertId,
int &  firstNonDirGraphVertID,
int &  nExtraDirichlet,
BottomUpSubStructuredGraphSharedPtr &  bottomUpGraph,
set< int > &  extraDirVerts,
const bool  checkIfSystemSingular = false,
int  mdswitch = 1,
bool  doInteriorMap = false 
)

Construct optimal ordering a two-dimensional expansion given a vector of boundary condition information

The only unique identifiers of the vertices and edges of the mesh are the vertex id and the mesh id (stored in their corresponding Geometry object). However, setting up a global numbering based on these id's will not lead to a suitable or optimal numbering. Mainly because:

The vertices and egdes therefore need to be rearranged which is perofrmed in in the following way: The vertices and edges of the mesh are considered as vertices of a graph (in a computer science terminology, equivalently, they can also be considered as boundary degrees of freedom, whereby all boundary modes of a single edge are considered as a single DOF). We then will use different algorithms to reorder the graph-vertices.

In the following we use a boost graph object to store this graph the first template parameter (=OutEdgeList) is chosen to be of type std::set. Similarly we also use a std::set to hold the adjacency information. A similar edge might exist multiple times and so to prevent the definition of parallel edges, we use std::set (=boost::setS) rather than std::vector (=boost::vecS).

Two different containers are used to store the graph vertex id's of the different mesh vertices and edges. They are implemented as a STL map such that the graph vertex id can later be retrieved by the unique mesh vertex or edge id's which serve as the key of the map.

Therefore, the algorithm proceeds as follows:

STEP 1: Order the Dirichlet vertices and edges first

STEP 1.4: Check for singular system and add pinning Dirichlet vertex

STEP 1.5: Exchange Dirichlet mesh vertices between processes and check for singular problems.

STEP 2: Now order all other vertices and edges in the graph and create a temporary numbering of domain-interior vertices and edges.

STEP 3: Reorder graph for optimisation.

STEP 4: Fill the #vertReorderedGraphVertId and #edgeReorderedGraphVertId with the optimal ordering from boost.

Definition at line 577 of file AssemblyMapCG2D.cpp.

References Nektar::StdRegions::StdExpansion::as(), ASSERTL0, ASSERTL1, Nektar::MultiRegions::CuthillMckeeReordering(), Nektar::MultiRegions::eDirectFullMatrix, Nektar::MultiRegions::eDirectMultiLevelStaticCond, Nektar::MultiRegions::eDirectStaticCond, Nektar::SpatialDomains::eDirichlet, Nektar::MultiRegions::eIterativeFull, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::MultiRegions::eIterativeStaticCond, Nektar::SpatialDomains::eNeumann, Nektar::MultiRegions::ePETScFullMatrix, Nektar::MultiRegions::ePETScMultiLevelStaticCond, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::eXxtFullMatrix, Nektar::MultiRegions::eXxtMultiLevelStaticCond, Nektar::MultiRegions::eXxtStaticCond, Gs::Gather(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::LocalRegions::Expansion1D::GetGeom1D(), Nektar::MultiRegions::ExpList::GetOffset_Elmt_Id(), Nektar::MultiRegions::GlobalSysSolnTypeMap, Gs::gs_add, Vmath::Imax(), Gs::Init(), Nektar::iterator, Nektar::MultiRegions::AssemblyMap::m_comm, Nektar::MultiRegions::AssemblyMap::m_lowestStaticCondLevel, Nektar::MultiRegions::AssemblyMap::m_numLocalBndCoeffs, Nektar::MultiRegions::AssemblyMapCG::m_numNonDirEdges, Nektar::MultiRegions::AssemblyMapCG::m_numNonDirVertexModes, Nektar::MultiRegions::AssemblyMap::m_session, Nektar::MultiRegions::AssemblyMap::m_solnType, Nektar::MultiRegions::MultiLevelBisectionReordering(), Nektar::MultiRegions::NoReordering(), Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceMin, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

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

{
int i,j,k,l;
int cnt = 0;
int meshVertId, meshVertId2;
int meshEdgeId;
int graphVertId = 0;
StdRegions::StdExpansion2DSharedPtr locExpansion;
LocalRegions::SegExpSharedPtr bndSegExp;
MultiRegions::ExpList0DSharedPtr bndVertExp;
const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
map<int,int>::iterator mapIt;
map<int,int>::const_iterator mapConstIt;
bool systemSingular = true;
LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
PeriodicMap::const_iterator pIt;
/**
* STEP 1: Order the Dirichlet vertices and edges first
*/
for(i = 0; i < bndCondExp.num_elements(); i++)
{
// Check to see if any value on edge has Dirichlet value.
cnt = 0;
for(k = 0; k < bndConditions.num_elements(); ++k)
{
if (bndConditions[k][i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
{
cnt++;
}
if (bndConditions[k][i]->GetBoundaryConditionType() !=
SpatialDomains::eNeumann)
{
systemSingular = false;
}
}
// If all boundaries are Dirichlet take out of mask
if(cnt == bndConditions.num_elements())
{
for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
{
bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
ReorderedGraphVertId[1][meshEdgeId] = graphVertId++;
for(k = 0; k < 2; k++)
{
meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
if(ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] =
graphVertId++;
}
}
}
}
}
/**
* STEP 1.4: Check for singular system and add pinning Dirichlet
* vertex
*/
// Check between processes if the whole system is singular
int n = m_comm->GetSize();
int p = m_comm->GetRank();
int s = systemSingular ? 1 : 0;
vCommRow->AllReduce(s, LibUtilities::ReduceMin);
systemSingular = (s == 1 ? true : false);
// Count the number of boundary regions on each process
Array<OneD, int> bccounts(n, 0);
bccounts[p] = bndCondExp.num_elements();
vCommRow->AllReduce(bccounts, LibUtilities::ReduceSum);
// Find the process rank with the maximum number of boundary regions
int maxBCIdx = Vmath::Imax(n, bccounts, 1);
// If the system is singular, the process with the maximum number of
// BCs will set a Dirichlet vertex to make system non-singular.
// Note: we find the process with maximum boundary regions to ensure
// we do not try to set a Dirichlet vertex on a partition with no
// intersection with the boundary.
meshVertId = 0;
if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p)
{
if (m_session->DefinesParameter("SingularElement"))
{
int s_eid;
m_session->LoadParameter("SingularElement", s_eid);
ASSERTL1(s_eid < locExpVector.size(),
"SingularElement Parameter is too large");
meshVertId = locExpVector[s_eid]
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetVid(0);
}
else if (m_session->DefinesParameter("SingularVertex"))
{
m_session->LoadParameter("SingularVertex", meshVertId);
}
else if (bndCondExp.num_elements() == 0)
{
// All boundaries are periodic.
meshVertId = locExpVector[0]
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetVid(0);
}
else
{
bndSegExp = bndCondExp[bndCondExp.num_elements()-1]
->GetExp(0)->as<LocalRegions::SegExp>();
meshVertId = bndSegExp->GetGeom1D()->GetVid(0);
}
if (ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] = graphVertId++;
}
}
vCommRow->AllReduce(meshVertId, LibUtilities::ReduceSum);
// When running in parallel, we need to ensure that the singular
// mesh vertex is communicated to any periodic vertices, otherwise
// the system may diverge.
if(systemSingular == true && checkIfSystemSingular)
{
// Firstly, we check that no other processors have this
// vertex. If they do, then we mark the vertex as also being
// Dirichlet.
if (maxBCIdx != p)
{
if (Dofs[0].count(meshVertId) > 0)
{
if (ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] = graphVertId++;
}
}
}
// In the case that meshVertId is periodic with other vertices,
// this process and all other processes need to make sure that
// the periodic vertices are also marked as Dirichlet.
int gId;
// At least one process (maxBCidx) will have already associated
// a graphVertId with meshVertId. Others won't even have any of
// the vertices. The logic below is designed to handle both
// cases.
if (ReorderedGraphVertId[0].count(meshVertId) == 0)
{
gId = -1;
}
else
{
gId = ReorderedGraphVertId[0][meshVertId];
}
for (pIt = periodicVerts.begin();
pIt != periodicVerts.end(); ++pIt)
{
// Either the vertex is local to this processor (in which
// case it will be in the pIt->first position) or else
// meshVertId might be contained within another processor's
// vertex list. The if statement below covers both cases. If
// we find it, set as Dirichlet with the vertex id gId.
if (pIt->first == meshVertId)
{
ReorderedGraphVertId[0][meshVertId] =
gId < 0 ? graphVertId++ : gId;
for (i = 0; i < pIt->second.size(); ++i)
{
if (pIt->second[i].isLocal)
{
ReorderedGraphVertId[0][pIt->second[i].id] =
gId;
}
}
}
else
{
bool found = false;
for (i = 0; i < pIt->second.size(); ++i)
{
if (pIt->second[i].id == meshVertId)
{
found = true;
break;
}
}
if (found)
{
ReorderedGraphVertId[0][pIt->first] =
gId < 0 ? graphVertId++ : gId;
for (i = 0; i < pIt->second.size(); ++i)
{
if (pIt->second[i].isLocal)
{
ReorderedGraphVertId[0][
pIt->second[i].id] = gId;
}
}
}
}
}
}
/**
* STEP 1.5: Exchange Dirichlet mesh vertices between processes and
* check for singular problems.
*/
// Collate information on Dirichlet vertices from all processes
Array<OneD, int> counts (n, 0);
Array<OneD, int> offsets(n, 0);
counts[p] = ReorderedGraphVertId[0].size();
vCommRow->AllReduce(counts, LibUtilities::ReduceSum);
for (i = 1; i < n; ++i)
{
offsets[i] = offsets[i-1] + counts[i-1];
}
int nTot = Vmath::Vsum(n,counts,1);
Array<OneD, int> vertexlist(nTot, 0);
std::map<int, int>::iterator it;
for (it = ReorderedGraphVertId[0].begin(), i = 0;
it != ReorderedGraphVertId[0].end();
++it, ++i)
{
vertexlist[offsets[p] + i] = it->first;
}
vCommRow->AllReduce(vertexlist, LibUtilities::ReduceSum);
map<int, int> extraDirVertIds;
// Ensure Dirchlet vertices are consistently recorded between
// processes (e.g. Dirichlet region meets Neumann region across a
// partition boundary requires vertex on partition to be Dirichlet).
for (i = 0; i < n; ++i)
{
if (i == p)
{
continue;
}
for(j = 0; j < locExpVector.size(); j++)
{
locExpansion = boost::dynamic_pointer_cast<
StdRegions::StdExpansion2D>(
locExpVector[locExp.GetOffset_Elmt_Id(j)]);
for(k = 0; k < locExpansion->GetNverts(); k++)
{
meshVertId = locExpansion
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetVid(k);
if (ReorderedGraphVertId[0].count(meshVertId) != 0)
{
continue;
}
for (l = 0; l < counts[i]; ++l)
{
if (vertexlist[offsets[i]+l] == meshVertId)
{
extraDirVertIds[meshVertId] = i;
ReorderedGraphVertId[0][meshVertId] =
graphVertId++;
nExtraDirichlet++;
}
}
}
}
}
for (i = 0; i < n; ++i)
{
counts [i] = 0;
offsets[i] = 0;
}
counts[p] = extraDirVertIds.size();
vCommRow->AllReduce(counts, LibUtilities::ReduceSum);
nTot = Vmath::Vsum(n, counts, 1);
offsets[0] = 0;
for (i = 1; i < n; ++i)
{
offsets[i] = offsets[i-1] + counts[i-1];
}
Array<OneD, int> vertids (nTot, 0);
Array<OneD, int> vertprocs(nTot, 0);
for (it = extraDirVertIds.begin(), i = 0;
it != extraDirVertIds.end(); ++it, ++i)
{
vertids [offsets[p]+i] = it->first;
vertprocs[offsets[p]+i] = it->second;
}
vCommRow->AllReduce(vertids, LibUtilities::ReduceSum);
vCommRow->AllReduce(vertprocs, LibUtilities::ReduceSum);
for (i = 0; i < nTot; ++i)
{
if (m_comm->GetRank() != vertprocs[i])
{
continue;
}
extraDirVerts.insert(vertids[i]);
}
firstNonDirGraphVertId = graphVertId;
/**
* STEP 2: Now order all other vertices and edges in the graph and
* create a temporary numbering of domain-interior vertices and
* edges.
*/
typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
BoostGraph boostGraphObj;
int tempGraphVertId = 0;
int nVerts;
int vertCnt;
int edgeCnt;
int localOffset=0;
int nTotalVerts=0;
map<int, int> vertTempGraphVertId;
map<int, int> edgeTempGraphVertId;
map<int, int> intTempGraphVertId;
Array<OneD, int> localVerts;
Array<OneD, int> localEdges;
Array<OneD, int> localinterior;
m_numLocalBndCoeffs = 0;
/// - Local periodic vertices.
for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
{
meshVertId = pIt->first;
// This periodic vertex is already Dirichlet.
if (ReorderedGraphVertId[0].count(pIt->first) != 0)
{
for (i = 0; i < pIt->second.size(); ++i)
{
meshVertId2 = pIt->second[i].id;
if (ReorderedGraphVertId[0].count(meshVertId2) == 0 &&
pIt->second[i].isLocal)
{
ReorderedGraphVertId[0][meshVertId2] =
ReorderedGraphVertId[0][meshVertId];
}
}
continue;
}
// One of the attached vertices is Dirichlet.
bool isDirichlet = false;
for (i = 0; i < pIt->second.size(); ++i)
{
if (!pIt->second[i].isLocal)
{
continue;
}
meshVertId2 = pIt->second[i].id;
if (ReorderedGraphVertId[0].count(meshVertId2) > 0)
{
isDirichlet = true;
break;
}
}
if (isDirichlet)
{
ReorderedGraphVertId[0][meshVertId] =
ReorderedGraphVertId[0][pIt->second[i].id];
for (j = 0; j < pIt->second.size(); ++j)
{
meshVertId2 = pIt->second[i].id;
if (j == i || !pIt->second[j].isLocal ||
ReorderedGraphVertId[0].count(meshVertId2) > 0)
{
continue;
}
ReorderedGraphVertId[0][meshVertId2] =
ReorderedGraphVertId[0][pIt->second[i].id];
}
continue;
}
// Otherwise, see if a vertex ID has already been set.
for (i = 0; i < pIt->second.size(); ++i)
{
if (!pIt->second[i].isLocal)
{
continue;
}
if (vertTempGraphVertId.count(pIt->second[i].id) > 0)
{
break;
}
}
if (i == pIt->second.size())
{
vertTempGraphVertId[meshVertId] = tempGraphVertId++;
}
else
{
vertTempGraphVertId[meshVertId] = vertTempGraphVertId[pIt->second[i].id];
}
}
/// - Local periodic edges. !! SJS: This will need pushing
/// - further down for block precodnitioner
for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
{
if (!pIt->second[0].isLocal)
{
// The edge mapped to is on another process.
meshEdgeId = pIt->first;
ASSERTL0(ReorderedGraphVertId[1].count(meshEdgeId) == 0,
"This periodic boundary edge has been specified before");
edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
}
else if (pIt->first < pIt->second[0].id)
{
ASSERTL0(ReorderedGraphVertId[1].count(pIt->first) == 0,
"This periodic boundary edge has been specified before");
ASSERTL0(ReorderedGraphVertId[1].count(pIt->second[0].id) == 0,
"This periodic boundary edge has been specified before");
edgeTempGraphVertId[pIt->first] = tempGraphVertId;
edgeTempGraphVertId[pIt->second[0].id] = tempGraphVertId++;
}
}
/// - All other vertices and edges
int elmtid;
for(i = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
locExpVector[elmtid])))
{
m_numLocalBndCoeffs += locExpansion->NumBndryCoeffs();
nTotalVerts += locExpansion->GetNverts();
}
}
// Store the temporary graph vertex
// id's of all element edges and
// vertices in these 2 arrays below
localVerts = Array<OneD, int>(nTotalVerts,-1);
localEdges = Array<OneD, int>(nTotalVerts,-1);
for(i = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = locExpVector[elmtid]
->as<StdRegions::StdExpansion2D>()))
{
vertCnt = 0;
nVerts = locExpansion->GetNverts();
for(j = 0; j < nVerts; ++j)
{
meshVertId = locExpansion
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetVid(j);
if(ReorderedGraphVertId[0].count(meshVertId) == 0)
{
// non-periodic & non-Dirichlet vertex
if(vertTempGraphVertId.count(meshVertId) == 0)
{
boost::add_vertex(boostGraphObj);
vertTempGraphVertId[meshVertId] = tempGraphVertId++;
}
localVerts[localOffset + vertCnt++] = vertTempGraphVertId[meshVertId];
}
}
}
localOffset+=nVerts;
}
m_numNonDirVertexModes=tempGraphVertId;
localOffset=0;
for(i = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = locExpVector[elmtid]
->as<StdRegions::StdExpansion2D>()))
{
edgeCnt = 0;
nVerts = locExpansion->GetNverts();
for(j = 0; j < nVerts; ++j)
{
meshEdgeId = locExpansion
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetEid(j);
if(ReorderedGraphVertId[1].count(meshEdgeId) == 0)
{
// non-periodic & non-Dirichlet edge
if(edgeTempGraphVertId.count(meshEdgeId) == 0)
{
boost::add_vertex(boostGraphObj);
edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
}
localEdges[localOffset + edgeCnt++] = edgeTempGraphVertId[meshEdgeId];
}
}
}
localOffset+=nVerts;
}
// Number of dirichlet edges
m_numNonDirEdges = edgeTempGraphVertId.size();
if(doInteriorMap)
{
for(i = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
locExpVector[elmtid])))
{
boost::add_vertex(boostGraphObj);
intTempGraphVertId[elmtid] = tempGraphVertId++;
}
}
}
localOffset=0;
for(i = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
locExpVector[elmtid])))
{
nVerts = locExpansion->GetNverts();
// Now loop over all local edges and vertices of this
// element and define that all other edges and vertices of
// this element are adjacent to them.
for(j = 0; j < nVerts; j++)
{
if(localVerts[j+localOffset]==-1)
{
break;
}
for(k = 0; k < nVerts; k++)
{
if(localVerts[k+localOffset]==-1)
{
break;
}
if(k!=j)
{
boost::add_edge( (size_t) localVerts[j+localOffset],
(size_t) localVerts[k+localOffset],boostGraphObj);
}
}
for(k = 0; k < nVerts; k++)
{
if(localEdges[k+localOffset]==-1)
{
break;
}
boost::add_edge( (size_t) localVerts[j+localOffset],
(size_t) localEdges[k+localOffset],boostGraphObj);
}
if(doInteriorMap)
{
boost::add_edge( (size_t) localVerts[j+localOffset],
(size_t) intTempGraphVertId[elmtid],boostGraphObj);
}
}
for(j = 0; j < nVerts; j++)
{
if(localEdges[j+localOffset]==-1)
{
break;
}
for(k = 0; k < nVerts; k++)
{
if(localEdges[k+localOffset]==-1)
{
break;
}
if(k!=j)
{
boost::add_edge( (size_t) localEdges[j+localOffset],
(size_t) localEdges[k+localOffset],boostGraphObj);
}
}
for(k = 0; k < nVerts; k++)
{
if(localVerts[k+localOffset]==-1)
{
break;
}
boost::add_edge( (size_t) localEdges[j+localOffset],
(size_t) localVerts[k+localOffset],boostGraphObj);
}
if(doInteriorMap)
{
boost::add_edge( (size_t) localEdges[j+localOffset],
(size_t) intTempGraphVertId[elmtid], boostGraphObj);
}
}
if(doInteriorMap)
{
for(j = 0; j < nVerts; j++)
{
if(localVerts[j+localOffset]==-1)
{
break;
}
boost::add_edge( (size_t) intTempGraphVertId[elmtid],
(size_t) localVerts[j+localOffset], boostGraphObj);
}
for(j = 0; j < nVerts; j++)
{
if(localEdges[j+localOffset]==-1)
{
break;
}
boost::add_edge( (size_t) intTempGraphVertId[elmtid],
(size_t) localEdges[j+localOffset], boostGraphObj);
}
}
}
else
{
ASSERTL0(false,"dynamic cast to a local 2D expansion failed");
}
localOffset+=nVerts;
}
// Container to store vertices of the graph which correspond to
// degrees of freedom along the boundary.
set<int> partVerts;
if (m_solnType == eIterativeMultiLevelStaticCond)
{
vector<long> procVerts, procEdges;
set <int> foundVerts, foundEdges;
// Loop over element and construct the procVerts and procEdges
// vectors, which store the geometry IDs of mesh vertices and
// edges respectively which are local to this process.
for(i = cnt = 0; i < locExpVector.size(); ++i)
{
elmtid = locExp.GetOffset_Elmt_Id(i);
if((locExpansion = boost::dynamic_pointer_cast<
StdRegions::StdExpansion2D>(locExpVector[elmtid])))
{
for (j = 0; j < locExpansion->GetNverts(); ++j, ++cnt)
{
int vid = locExpansion
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetVid(j)+1;
int eid = locExpansion
->as<LocalRegions::Expansion2D>()
->GetGeom2D()->GetEid(j)+1;
if (foundVerts.count(vid) == 0)
{
procVerts.push_back(vid);
foundVerts.insert(vid);
}
if (foundEdges.count(eid) == 0)
{
procEdges.push_back(eid);
foundEdges.insert(eid);
}
}
}
else
{
ASSERTL0(false,
"dynamic cast to a local 2D expansion failed");
}
}
int unique_verts = foundVerts.size();
int unique_edges = foundEdges.size();
// Now construct temporary GS objects. These will be used to
// populate the arrays tmp3 and tmp4 with the multiplicity of
// the vertices and edges respectively to identify those
// vertices and edges which are located on partition boundary.
Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
Gs::gs_data *tmp1 = Gs::Init(vertArray, m_comm);
Gs::gs_data *tmp2 = Gs::Init(edgeArray, m_comm);
Array<OneD, NekDouble> tmp3(unique_verts, 1.0);
Array<OneD, NekDouble> tmp4(unique_edges, 1.0);
Gs::Gather(tmp3, Gs::gs_add, tmp1);
Gs::Gather(tmp4, Gs::gs_add, tmp2);
// Finally, fill the partVerts set with all non-Dirichlet
// vertices which lie on a partition boundary.
for (i = 0; i < unique_verts; ++i)
{
if (tmp3[i] > 1.0)
{
if (ReorderedGraphVertId[0].count(procVerts[i]-1) == 0)
{
partVerts.insert(vertTempGraphVertId[procVerts[i]-1]);
}
}
}
for (i = 0; i < unique_edges; ++i)
{
if (tmp4[i] > 1.0)
{
if (ReorderedGraphVertId[1].count(procEdges[i]-1) == 0)
{
partVerts.insert(edgeTempGraphVertId[procEdges[i]-1]);
}
}
}
}
/**
* STEP 3: Reorder graph for optimisation.
*/
int nGraphVerts = tempGraphVertId;
Array<OneD, int> perm(nGraphVerts);
Array<OneD, int> iperm(nGraphVerts);
if(nGraphVerts)
{
switch(m_solnType)
{
case eDirectFullMatrix:
case eIterativeFull:
case eIterativeStaticCond:
case eXxtFullMatrix:
case eXxtStaticCond:
case ePETScFullMatrix:
case ePETScStaticCond:
{
NoReordering(boostGraphObj,perm,iperm);
}
break;
case eDirectStaticCond:
{
CuthillMckeeReordering(boostGraphObj,perm,iperm);
}
break;
case eDirectMultiLevelStaticCond:
case eIterativeMultiLevelStaticCond:
case eXxtMultiLevelStaticCond:
case ePETScMultiLevelStaticCond:
{
MultiLevelBisectionReordering(boostGraphObj,perm,iperm,bottomUpGraph,partVerts,mdswitch);
}
break;
default:
{
ASSERTL0(false,"Unrecognised solution type: " + std::string(MultiRegions::GlobalSysSolnTypeMap[m_solnType]));
}
}
}
// For parallel multi-level static condensation determine the lowest
// static condensation level amongst processors.
if (m_solnType == eIterativeMultiLevelStaticCond)
{
m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
vCommRow->AllReduce(m_lowestStaticCondLevel,
LibUtilities::ReduceMax);
}
else
{
m_lowestStaticCondLevel = 0;
}
/**
* STEP 4: Fill the #vertReorderedGraphVertId and
* #edgeReorderedGraphVertId with the optimal ordering from boost.
*/
for(mapIt = vertTempGraphVertId.begin(); mapIt != vertTempGraphVertId.end(); mapIt++)
{
ReorderedGraphVertId[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
}
for(mapIt = edgeTempGraphVertId.begin(); mapIt != edgeTempGraphVertId.end(); mapIt++)
{
ReorderedGraphVertId[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
}
if(doInteriorMap)
{
for(mapIt = intTempGraphVertId.begin(); mapIt != intTempGraphVertId.end(); mapIt++)
{
ReorderedGraphVertId[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
}
}
return nGraphVerts;
}
Generated on Sun Mar 15 2015 00:11:57 for Nektar++ by   doxygen 1.8.1.2