Constructs mappings for the C0 scalar continuous Galerkin formulation. More...
#include <AssemblyMapCG.h>
Inheritance diagram for Nektar::MultiRegions::AssemblyMapCG:
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... | |
| ~AssemblyMapCG () override | |
| 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... | |
| std::string | GetVariable () |
| Retrieves the variable string. 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... | |
| std::string | GetPreconType () const |
| NekDouble | GetIterativeTolerance () const |
| bool | IsAbsoluteTolerance () 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... | |
| int | v_GetLocalToGlobalMap (const int i) const override |
| int | v_GetGlobalToUniversalMap (const int i) const override |
| int | v_GetGlobalToUniversalMapUnique (const int i) const override |
| const Array< OneD, const int > & | v_GetLocalToGlobalMap () override |
| const Array< OneD, const int > & | v_GetGlobalToUniversalMap () override |
| const Array< OneD, const int > & | v_GetGlobalToUniversalMapUnique () override |
| NekDouble | v_GetLocalToGlobalSign (const int i) const override |
| const Array< OneD, NekDouble > & | v_GetLocalToGlobalSign () const override |
| void | v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const override |
| void | v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm) const override |
| void | v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const override |
| void | v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const override |
| void | v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const override |
| void | v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const override |
| void | v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const override |
| void | v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const override |
| void | v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const override |
| int | v_GetFullSystemBandWidth () const override |
| int | v_GetNumNonDirVertexModes () const override |
| int | v_GetNumNonDirEdgeModes () const override |
| int | v_GetNumNonDirFaceModes () const override |
| int | v_GetNumDirEdges () const override |
| int | v_GetNumDirFaces () const override |
| int | v_GetNumNonDirEdges () const override |
| int | v_GetNumNonDirFaces () const override |
| const Array< OneD, const int > & | v_GetExtraDirEdges () override |
| 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) |
| virtual int | v_GetLocalToGlobalMap (const int i) const |
| virtual int | v_GetGlobalToUniversalMap (const int i) const |
| virtual int | v_GetGlobalToUniversalMapUnique (const int i) const |
| virtual const Array< OneD, const int > & | v_GetLocalToGlobalMap () |
| virtual const Array< OneD, const int > & | v_GetGlobalToUniversalMap () |
| virtual const Array< OneD, const int > & | v_GetGlobalToUniversalMapUnique () |
| virtual NekDouble | v_GetLocalToGlobalSign (const int i) const |
| virtual const Array< OneD, NekDouble > & | v_GetLocalToGlobalSign () const |
| virtual void | v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const |
| virtual void | v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, bool useComm) const |
| virtual void | v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const |
| virtual void | v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const |
| virtual void | v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const |
| virtual void | v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const |
| virtual void | v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const |
| virtual void | v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const |
| virtual void | v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const |
| virtual int | v_GetFullSystemBandWidth () const |
| virtual int | v_GetNumNonDirVertexModes () const |
| virtual int | v_GetNumNonDirEdgeModes () const |
| virtual int | v_GetNumNonDirFaceModes () const |
| virtual int | v_GetNumDirEdges () const |
| virtual int | v_GetNumDirFaces () const |
| virtual int | v_GetNumNonDirEdges () const |
| virtual int | v_GetNumNonDirFaces () const |
| virtual const Array< OneD, const int > & | v_GetExtraDirEdges () |
| virtual std::shared_ptr< AssemblyMap > | v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType) |
| Generate a linear space mapping from existing mapping. More... | |
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... | |
| std::string | m_variable |
| Variable string identifier. 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... | |
| std::string | m_preconType |
| Type type of preconditioner to use in iterative solver. More... | |
| NekDouble | m_iterativeTolerance |
| Tolerance for iterative solver. More... | |
| bool | m_isAbsoluteTolerance |
| 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:
- Direct full matrix
- Direct static condensation
- Direct multi-level static condensation In the latter case, mappings are created recursively for the different levels of static condensation.
These mappings are used by GlobalLinSys to generate the global system.
Definition at line 65 of file AssemblyMapCG.h.
Member Typedef Documentation
◆ BndCond
|
private |
Definition at line 68 of file AssemblyMapCG.h.
◆ BndCondExp
|
private |
Definition at line 67 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 69 of file AssemblyMapCG.cpp.
72 : AssemblyMap(pSession, comm, variable)
73{
75}
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 1319 of file AssemblyMapCG.cpp.
1326 : AssemblyMap(pSession, locExp.GetComm(), variable)
1327{
1328 int i, j, k;
1330 int cnt = 0;
1331 int meshVertId, meshEdgeId, meshEdgeId2, meshFaceId, meshFaceId2;
1332 int globalId;
1333 int nEdgeInteriorCoeffs;
1334 int firstNonDirGraphVertId;
1336 LocalRegions::ExpansionSharedPtr exp, bndExp;
1337 StdRegions::Orientation edgeOrient;
1338 StdRegions::Orientation faceOrient;
1339 Array<OneD, unsigned int> edgeInteriorMap;
1340 Array<OneD, int> edgeInteriorSign;
1341 Array<OneD, unsigned int> faceInteriorMap;
1342 Array<OneD, int> faceInteriorSign;
1343
1345
1347
1349
1350 // Stores vertex, edge and face reordered vertices.
1351 DofGraph graph(3);
1352 DofGraph dofs(3);
1353 vector<map<int, int>> faceModes(2);
1354 map<int, LibUtilities::ShapeType> faceType;
1355
1356 set<int> extraDirVerts, extraDirEdges;
1357 BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
1358
1359 // Construct list of number of degrees of freedom for each vertex,
1360 // edge and face.
1361 for (i = 0; i < locExpVector.size(); ++i)
1362 {
1363 exp = locExpVector[i];
1364
1365 for (j = 0; j < exp->GetNverts(); ++j)
1366 {
1367 dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1368 }
1369
1370 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1371 {
1372 int nEdgeInt;
1373 if (exp->GetGeom()->GetNumFaces())
1374 {
1375 nEdgeInt =
1376 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
1377 }
1378 else
1379 {
1380 nEdgeInt = exp->GetTraceNcoeffs(j) - 2;
1381 }
1382
1383 if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1384 {
1385 if (dofs[1][exp->GetGeom()->GetEid(j)] != nEdgeInt)
1386 {
1387 ASSERTL0(
1388 (exp->GetBasisType(0) == LibUtilities::eModified_A) ||
1389 (exp->GetBasisType(1) ==
1390 LibUtilities::eModified_B) ||
1391 (exp->GetBasisType(2) ==
1392 LibUtilities::eModified_C) ||
1393 (exp->GetBasisType(2) ==
1394 LibUtilities::eModifiedPyr_C),
1395 "CG with variable order only available with "
1396 "modal expansion");
1397 }
1398 dofs[1][exp->GetGeom()->GetEid(j)] =
1399 min(dofs[1][exp->GetGeom()->GetEid(j)], nEdgeInt);
1400 }
1401 else
1402 {
1403 dofs[1][exp->GetGeom()->GetEid(j)] = nEdgeInt;
1404 }
1405 }
1406
1407 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1408 {
1409 faceOrient = exp->GetGeom()->GetForient(j);
1410 meshFaceId = exp->GetGeom()->GetFid(j);
1411 exp->GetTraceNumModes(j, numModes0, numModes1, faceOrient);
1412
1413 if (faceModes[0].count(meshFaceId) > 0)
1414 {
1415 faceModes[0][meshFaceId] =
1416 min(faceModes[0][meshFaceId], numModes0);
1417
1418 faceModes[1][meshFaceId] =
1419 min(faceModes[1][meshFaceId], numModes1);
1420 }
1421 else
1422 {
1423 faceModes[0][meshFaceId] = numModes0;
1424 faceModes[1][meshFaceId] = numModes1;
1425
1426 // Get shape of this face
1427 SpatialDomains::Geometry3DSharedPtr geom;
1428 geom = std::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
1429 exp->GetGeom());
1430 faceType[meshFaceId] = geom->GetFace(j)->GetShapeType();
1431 }
1432 }
1433 }
1434
1435 // Add non-local periodic dofs to the map
1436 for (auto &pIt : periodicEdges)
1437 {
1438 for (i = 0; i < pIt.second.size(); ++i)
1439 {
1440 meshEdgeId2 = pIt.second[i].id;
1441 if (dofs[1].count(meshEdgeId2) == 0)
1442 {
1443 dofs[1][meshEdgeId2] = 1e6;
1444 }
1445 }
1446 }
1447 for (auto &pIt : periodicFaces)
1448 {
1449 for (i = 0; i < pIt.second.size(); ++i)
1450 {
1451 meshFaceId2 = pIt.second[i].id;
1452 if (faceModes[0].count(meshFaceId2) == 0)
1453 {
1454 faceModes[0][meshFaceId2] = 1e6;
1455 faceModes[1][meshFaceId2] = 1e6;
1456 }
1457 }
1458 }
1459
1460 // Now use information from all partitions to determine the correct
1461 // size
1462
1463 // edges
1464 Array<OneD, long> edgeId(dofs[1].size());
1465 Array<OneD, NekDouble> edgeDof(dofs[1].size());
1466 i = 0;
1467 for (auto &dofIt : dofs[1])
1468 {
1469 edgeId[i] = dofIt.first + 1;
1470 edgeDof[i++] = (NekDouble)dofIt.second;
1471 }
1474 Gs::Finalise(tmp);
1475 for (i = 0; i < dofs[1].size(); i++)
1476 {
1477 dofs[1][edgeId[i] - 1] = (int)(edgeDof[i] + 0.5);
1478 }
1479 // Periodic edges
1480 for (auto &pIt : periodicEdges)
1481 {
1482 meshEdgeId = pIt.first;
1483 for (i = 0; i < pIt.second.size(); ++i)
1484 {
1485 meshEdgeId2 = pIt.second[i].id;
1486 if (dofs[1][meshEdgeId2] < dofs[1][meshEdgeId])
1487 {
1488 dofs[1][meshEdgeId] = dofs[1][meshEdgeId2];
1489 }
1490 }
1491 }
1492 // faces
1493 Array<OneD, long> faceId(faceModes[0].size());
1494 Array<OneD, NekDouble> faceP(faceModes[0].size());
1495 Array<OneD, NekDouble> faceQ(faceModes[0].size());
1496
1497 i = 0;
1498 for (auto dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin();
1499 dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
1500 {
1501 faceId[i] = dofIt->first + 1;
1502 faceP[i] = (NekDouble)dofIt->second;
1503 faceQ[i] = (NekDouble)dofIt2->second;
1504 }
1508 Gs::Finalise(tmp2);
1509 for (i = 0; i < faceModes[0].size(); i++)
1510 {
1511 faceModes[0][faceId[i] - 1] = (int)(faceP[i] + 0.5);
1512 faceModes[1][faceId[i] - 1] = (int)(faceQ[i] + 0.5);
1513 }
1514 // Periodic faces
1515 for (auto &pIt : periodicFaces)
1516 {
1517 meshFaceId = pIt.first;
1518 for (i = 0; i < pIt.second.size(); ++i)
1519 {
1520 meshFaceId2 = pIt.second[i].id;
1521 if (faceModes[0][meshFaceId2] < faceModes[0][meshFaceId])
1522 {
1523 faceModes[0][meshFaceId] = faceModes[0][meshFaceId2];
1524 }
1525 if (faceModes[1][meshFaceId2] < faceModes[1][meshFaceId])
1526 {
1527 faceModes[1][meshFaceId] = faceModes[1][meshFaceId2];
1528 }
1529 }
1530 }
1531 // Calculate number of dof in each face
1533 for (i = 0; i < faceModes[0].size(); i++)
1534 {
1535 P = faceModes[0][faceId[i] - 1];
1536 Q = faceModes[1][faceId[i] - 1];
1538 {
1539 // Quad face
1540 dofs[2][faceId[i] - 1] =
1543 }
1544 else
1545 {
1546 // Tri face
1547 dofs[2][faceId[i] - 1] =
1550 }
1551 }
1552
1553 Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1554
1555 // Note that nExtraDirichlet is not used in the logic below; it just
1556 // needs to be set so that the coupled solver in
1557 // IncNavierStokesSolver can work.
1558 int nExtraDirichlet;
1559 int mdswitch;
1561
1563 locExp, bndCondExp, bndCondVec, checkIfSystemSingular, periodicVerts,
1564 periodicEdges, periodicFaces, graph, bottomUpGraph, extraDirVerts,
1565 extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, mdswitch);
1566
1567 /*
1568 * Set up an array which contains the offset information of the
1569 * different graph vertices.
1570 *
1571 * This basically means to identify to how many global degrees of
1572 * freedom the individual graph vertices correspond. Obviously,
1573 * the graph vertices corresponding to the mesh-vertices account
1574 * for a single global DOF. However, the graph vertices
1575 * corresponding to the element edges correspond to N-2 global DOF
1576 * where N is equal to the number of boundary modes on this edge.
1577 */
1578 Array<OneD, int> graphVertOffset(
1579 graph[0].size() + graph[1].size() + graph[2].size() + 1, 0);
1580
1581 graphVertOffset[0] = 0;
1582
1583 for (i = 0; i < locExpVector.size(); ++i)
1584 {
1585 exp = locExpVector[i];
1586
1587 for (j = 0; j < exp->GetNverts(); ++j)
1588 {
1589 meshVertId = exp->GetGeom()->GetVid(j);
1590 graphVertOffset[graph[0][meshVertId] + 1] = 1;
1591 }
1592
1593 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1594 {
1595 if (exp->GetGeom()->GetNumFaces()) // 3D version
1596 {
1597 nEdgeInteriorCoeffs =
1598 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
1599 }
1600 else
1601 {
1602 nEdgeInteriorCoeffs = exp->GetTraceNcoeffs(j) - 2;
1603 }
1604 meshEdgeId = exp->GetGeom()->GetEid(j);
1605 graphVertOffset[graph[1][meshEdgeId] + 1] = dofs[1][meshEdgeId];
1606
1607 // Need a sign vector for modal expansions if nEdgeCoeffs
1608 // >=3 (not 4 because of variable order case)
1609 if (nEdgeInteriorCoeffs &&
1610 (exp->GetBasisType(0) == LibUtilities::eModified_A))
1611 {
1613 }
1614 }
1615
1616 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1617 {
1618 meshFaceId = exp->GetGeom()->GetFid(j);
1619 graphVertOffset[graph[2][meshFaceId] + 1] = dofs[2][meshFaceId];
1620 }
1621 }
1622
1623 for (i = 1; i < graphVertOffset.size(); i++)
1624 {
1625 graphVertOffset[i] += graphVertOffset[i - 1];
1626 }
1627
1628 // Allocate the proper amount of space for the class-data
1629 m_numLocalCoeffs = numLocalCoeffs;
1630 m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1634 m_localToLocalIntMap =
1637 Array<OneD, int>(m_numLocalBndCondCoeffs, -1);
1638
1639 // If required, set up the sign-vector
1641 {
1643 m_localToGlobalBndSign =
1644 Array<OneD, NekDouble>(m_numLocalBndCoeffs, 1.0);
1646 Array<OneD, NekDouble>(m_numLocalBndCondCoeffs, 1.0);
1647 }
1648
1649 m_staticCondLevel = 0;
1650 m_numPatches = locExpVector.size();
1654 {
1655 m_numLocalBndCoeffsPerPatch[i] =
1656 (unsigned int)locExpVector[i]->NumBndryCoeffs();
1657 m_numLocalIntCoeffsPerPatch[i] =
1658 (unsigned int)locExpVector[i]->GetNcoeffs() -
1659 locExpVector[i]->NumBndryCoeffs();
1660 }
1661
1662 /**
1663 * STEP 6: Now, all ingredients are ready to set up the actual
1664 * local to global mapping.
1665 *
1666 * The remainder of the map consists of the element-interior
1667 * degrees of freedom. This leads to the block-diagonal submatrix
1668 * as each element-interior mode is globally orthogonal to modes
1669 * in all other elements.
1670 */
1671 cnt = 0;
1672
1673 // Loop over all the elements in the domain
1674 int cntbdry = 0;
1675 int cntint = 0;
1676 for (i = 0; i < locExpVector.size(); ++i)
1677 {
1678 exp = locExpVector[i];
1679 cnt = locExp.GetCoeff_Offset(i);
1680
1681 int nbdry = exp->NumBndryCoeffs();
1682 int nint = exp->GetNcoeffs() - nbdry;
1683
1684 Array<OneD, unsigned int> bmap(nbdry);
1685 Array<OneD, unsigned int> imap(nint);
1686
1687 exp->GetBoundaryMap(bmap);
1688 exp->GetInteriorMap(imap);
1689
1690 for (j = 0; j < nbdry; ++j)
1691 {
1692 m_localToLocalBndMap[cntbdry++] = cnt + bmap[j];
1693 }
1694
1695 for (j = 0; j < nint; ++j)
1696 {
1697 m_localToLocalIntMap[cntint++] = cnt + imap[j];
1698 }
1699
1700 for (j = 0; j < exp->GetNverts(); ++j)
1701 {
1702 meshVertId = exp->GetGeom()->GetVid(j);
1703
1704 // Set the global DOF for vertex j of element i
1705 m_localToGlobalMap[cnt + exp->GetVertexMap(j)] =
1706 graphVertOffset[graph[0][meshVertId]];
1707 }
1708
1709 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1710 {
1711 if (exp->GetGeom()->GetNumFaces())
1712 {
1713 nEdgeInteriorCoeffs =
1714 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
1715 }
1716 else
1717 {
1718 nEdgeInteriorCoeffs = exp->GetTraceNcoeffs(j) - 2;
1719 }
1720 edgeOrient = exp->GetGeom()->GetEorient(j);
1721 meshEdgeId = exp->GetGeom()->GetEid(j);
1722
1723 auto pIt = periodicEdges.find(meshEdgeId);
1724
1725 // See if this edge is periodic. If it is, then we map all
1726 // edges to the one with lowest ID, and align all
1727 // coefficients to this edge orientation.
1728 if (pIt != periodicEdges.end())
1729 {
1730 pair<int, StdRegions::Orientation> idOrient =
1731 DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
1732 pIt->second);
1733 edgeOrient = idOrient.second;
1734 }
1735
1736 if (exp->GetGeom()->GetNumFaces())
1737 {
1738 exp->as<LocalRegions::Expansion3D>()
1739 ->GetEdgeInteriorToElementMap(j, edgeInteriorMap,
1740 edgeInteriorSign, edgeOrient);
1741 }
1742 else
1743 {
1744 exp->GetTraceInteriorToElementMap(j, edgeInteriorMap,
1745 edgeInteriorSign, edgeOrient);
1746 }
1747
1748 // Set the global DOF's for the interior modes of edge j
1749 for (k = 0; k < dofs[1][meshEdgeId]; ++k)
1750 {
1751 m_localToGlobalMap[cnt + edgeInteriorMap[k]] =
1752 graphVertOffset[graph[1][meshEdgeId]] + k;
1753 }
1754 for (k = dofs[1][meshEdgeId]; k < nEdgeInteriorCoeffs; ++k)
1755 {
1756 m_localToGlobalMap[cnt + edgeInteriorMap[k]] = 0;
1757 }
1758
1759 // Fill the sign vector if required
1761 {
1762 for (k = 0; k < dofs[1][meshEdgeId]; ++k)
1763 {
1764 m_localToGlobalSign[cnt + edgeInteriorMap[k]] =
1765 (NekDouble)edgeInteriorSign[k];
1766 }
1767 for (k = dofs[1][meshEdgeId]; k < nEdgeInteriorCoeffs; ++k)
1768 {
1769 m_localToGlobalSign[cnt + edgeInteriorMap[k]] = 0.0;
1770 }
1771 }
1772 }
1773
1774 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1775 {
1776 faceOrient = exp->GetGeom()->GetForient(j);
1777 meshFaceId = exp->GetGeom()->GetFid(j);
1778
1779 auto pIt = periodicFaces.find(meshFaceId);
1780
1781 if (pIt != periodicFaces.end() &&
1782 meshFaceId == min(meshFaceId, pIt->second[0].id))
1783 {
1784 faceOrient = DeterminePeriodicFaceOrient(faceOrient,
1785 pIt->second[0].orient);
1786 }
1787
1788 exp->GetTraceInteriorToElementMap(j, faceInteriorMap,
1789 faceInteriorSign, faceOrient);
1790
1791 // Set the global DOF's for the interior modes of face j
1792 exp->GetTraceNumModes(j, numModes0, numModes1, faceOrient);
1793 switch (faceType[meshFaceId])
1794 {
1796 {
1797 int kLoc = 0;
1798 k = 0;
1800 {
1802 {
1804 (q < faceModes[1][meshFaceId]))
1805 {
1806 m_localToGlobalMap[cnt +
1807 faceInteriorMap[kLoc]] =
1808 graphVertOffset[graph[2][meshFaceId]] + k;
1810 {
1811 m_localToGlobalSign[cnt +
1812 faceInteriorMap[kLoc]] =
1813 (NekDouble)faceInteriorSign[kLoc];
1814 }
1815 k++;
1816 }
1817 else
1818 {
1819 m_localToGlobalMap[cnt +
1820 faceInteriorMap[kLoc]] = 0;
1822 {
1823 m_localToGlobalSign[cnt +
1824 faceInteriorMap[kLoc]] =
1825 0.0;
1826 }
1827 }
1828 kLoc++;
1829 }
1830 }
1831 }
1832 break;
1834 {
1835 int kLoc = 0;
1836 k = 0;
1838 {
1840 {
1843 {
1844 m_localToGlobalMap[cnt +
1845 faceInteriorMap[kLoc]] =
1846 graphVertOffset[graph[2][meshFaceId]] + k;
1848 {
1849 m_localToGlobalSign[cnt +
1850 faceInteriorMap[kLoc]] =
1851 (NekDouble)faceInteriorSign[kLoc];
1852 }
1853 k++;
1854 }
1855 else
1856 {
1857 m_localToGlobalMap[cnt +
1858 faceInteriorMap[kLoc]] = 0;
1860 {
1861 m_localToGlobalSign[cnt +
1862 faceInteriorMap[kLoc]] =
1863 0.0;
1864 }
1865 }
1866 kLoc++;
1867 }
1868 }
1869 }
1870 break;
1871 default:
1873 break;
1874 }
1875 }
1876 }
1877
1878 // Set up the mapping for the boundary conditions
1879 // Set up boundary mapping
1880 map<int, pair<int, int>> traceToElmtTraceMap;
1881 int id;
1882
1883 for (cnt = i = 0; i < locExpVector.size(); ++i)
1884 {
1885 exp = locExpVector[i];
1886
1887 for (j = 0; j < exp->GetNtraces(); ++j)
1888 {
1889 id = exp->GetGeom()->GetTid(j);
1890
1891 traceToElmtTraceMap[id] = pair<int, int>(i, j);
1892 }
1893 }
1894
1895 Array<OneD, unsigned int> maparray;
1896 Array<OneD, int> signarray;
1897 map<int, pair<int, NekDouble>> GloDirBndCoeffToLocalCoeff;
1898 set<int> CoeffOnDirTrace;
1899
1900 cnt = 0;
1901 int offset = 0;
1902 for (i = 0; i < bndCondExp.size(); i++)
1903 {
1904 set<int> foundExtraVerts, foundExtraEdges;
1905 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1906 {
1907 bndExp = bndCondExp[i]->GetExp(j);
1908 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1909
1910 int id = bndExp->GetGeom()->GetGlobalID();
1911
1913 "Failed to find trace id");
1914
1915 int eid = traceToElmtTraceMap[id].first;
1916 int tid = traceToElmtTraceMap[id].second;
1917
1918 exp = locExpVector[eid];
1919 int dim = exp->GetShapeDimension();
1920
1921 if (dim == 1)
1922 {
1923 m_bndCondCoeffsToLocalCoeffsMap[cnt] =
1924 locExp.GetCoeff_Offset(eid) + exp->GetVertexMap(tid);
1925 }
1926 else
1927 {
1928 if (dim == 2)
1929 {
1930 exp->GetTraceToElementMap(tid, maparray, signarray,
1931 exp->GetGeom()->GetEorient(tid),
1932 bndExp->GetBasisNumModes(0));
1933 }
1934 else if (dim == 3)
1935 {
1936 exp->GetTraceToElementMap(tid, maparray, signarray,
1937 exp->GetGeom()->GetForient(tid),
1938 bndExp->GetBasisNumModes(0),
1939 bndExp->GetBasisNumModes(1));
1940 }
1941
1942 for (k = 0; k < bndExp->GetNcoeffs(); k++)
1943 {
1944 m_bndCondCoeffsToLocalCoeffsMap[cnt + k] =
1945 locExp.GetCoeff_Offset(eid) + maparray[k];
1947 {
1948 m_bndCondCoeffsToLocalCoeffsSign[cnt + k] =
1949 signarray[k];
1950 }
1951 }
1952 }
1953
1954 // we now need some information to work out how to
1955 // handle vertices and edges that are only just
1956 // touching a dirichlet boundary (and not the
1957 // whole edge/face)
1958
1959 for (k = 0; k < bndExp->GetNcoeffs(); k++)
1960 {
1964
1966 {
1968 }
1969
1970 if (bndConditions[i]->GetBoundaryConditionType() ==
1972 {
1973 CoeffOnDirTrace.insert(locid);
1974
1975 // store the local id and sign from global id
1976 // back to local space;
1977 GloDirBndCoeffToLocalCoeff[gloid] =
1978 pair<int, NekDouble>(locid, sign);
1979 }
1980 }
1981 }
1982 offset += bndCondExp[i]->GetNcoeffs();
1983 }
1984
1986 m_numGlobalBndCoeffs = globalId;
1987
1988 // Set up a mapping list of Dirichlet Local Dofs that
1989 // arise due to one vertex or edge just touching a
1990 // Dirichlet boundary and need the value from another
1991 // local coeff that has been filled by the boundary
1992 // coeffs.
1993
1994 Array<OneD, NekDouble> gloParaDirBnd(m_numGlobalBndCoeffs, -1.0);
1995
1996 Array<OneD, unsigned int> bndmap;
1997 cnt = 0;
1998 for (i = 0; i < locExpVector.size(); ++i)
1999 {
2000 int gloid;
2001
2002 exp = locExpVector[i];
2003
2004 exp->GetBoundaryMap(bndmap);
2005
2006 for (j = 0; j < bndmap.size(); ++j)
2007 {
2008 k = cnt + bndmap[j];
2009
2010 if (CoeffOnDirTrace.count(k) == 0)
2011 {
2012 gloid = m_localToGlobalMap[k];
2013
2015 {
2016 if (GloDirBndCoeffToLocalCoeff.count(gloid))
2017 {
2018 int locid = GloDirBndCoeffToLocalCoeff[gloid].first;
2020
2022 {
2024 m_localToGlobalSign[k];
2025 }
2026
2028 // could make same `structure as extraDirDof
2029 m_copyLocalDirDofs.insert(DirDofs);
2030 }
2031 else // else could be on another parallel partition.
2032 {
2033 gloParaDirBnd[gloid] = gloid;
2034 }
2035 }
2036 }
2037 }
2038
2039 cnt += exp->GetNcoeffs();
2040 }
2041
2042 /*
2043 * The boundary condition mapping is generated from the same vertex
2044 * renumbering.
2045 */
2046 cnt = 0;
2048 {
2050 {
2051 m_localToGlobalMap[i] = globalId++;
2052 }
2053 else
2054 {
2056 {
2058 }
2060 }
2061 }
2062
2063 m_numGlobalCoeffs = globalId;
2064
2065 SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges,
2066 periodicFaces);
2067
2068 // Now that universal map is setup reset gloParaDirBnd to
2069 // 0 if no point communicated or universal value of not
2070 // equal to -1.0
2072 {
2073 int gloid = gloParaDirBnd[i];
2074 if (gloid == -1)
2075 {
2076 gloParaDirBnd[i] = 0.0;
2077 }
2078 else
2079 {
2080 gloParaDirBnd[i] = m_globalToUniversalMap[gloid];
2081 }
2082 }
2083
2084 // Use parallel boundary communication to set parallel
2085 // dirichlet values on all processors Needs to be after
2086 // SetupUuniversialC0ContMap
2088
2089 // copy global ids back to local values in partition to
2090 // initialise gs communicator.
2091 Array<OneD, long> paraDirBnd(m_numLocalCoeffs);
2092 for (i = 0; i < numLocalCoeffs; ++i)
2093 {
2094 paraDirBnd[i] = 0.0;
2095
2097
2099 {
2100 continue;
2101 }
2102
2103 paraDirBnd[i] = gloParaDirBnd[id];
2104
2105 if (gloParaDirBnd[id] > 0.0)
2106 {
2107 // gather any sign changes due to edge modes
2109 {
2111 {
2112 m_parallelDirBndSign.insert(i);
2113 }
2114 }
2115 }
2116 }
2117
2119
2120 // Set up the local to global map for the next level when using
2121 // multi-level static condensation
2126 nGraphVerts)
2127 {
2129 {
2130 Array<OneD, int> vwgts_perm(graph[0].size() + graph[1].size() +
2131 graph[2].size() -
2132 firstNonDirGraphVertId);
2133
2134 for (i = 0; i < locExpVector.size(); ++i)
2135 {
2136 exp = locExpVector[i];
2137
2138 for (j = 0; j < exp->GetNverts(); ++j)
2139 {
2140 meshVertId = exp->GetGeom()->GetVid(j);
2141
2142 if (graph[0][meshVertId] >= firstNonDirGraphVertId)
2143 {
2144 vwgts_perm[graph[0][meshVertId] -
2145 firstNonDirGraphVertId] =
2146 dofs[0][meshVertId];
2147 }
2148 }
2149
2150 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
2151 {
2152 meshEdgeId = exp->GetGeom()->GetEid(j);
2153
2154 if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
2155 {
2156 vwgts_perm[graph[1][meshEdgeId] -
2157 firstNonDirGraphVertId] =
2158 dofs[1][meshEdgeId];
2159 }
2160 }
2161
2162 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
2163 {
2164 meshFaceId = exp->GetGeom()->GetFid(j);
2165
2166 if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
2167 {
2168 vwgts_perm[graph[2][meshFaceId] -
2169 firstNonDirGraphVertId] =
2170 dofs[2][meshFaceId];
2171 }
2172 }
2173 }
2174
2175 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
2176 m_nextLevelLocalToGlobalMap =
2178 bottomUpGraph);
2179 }
2180 }
2181
2183
2184 // Add up hash values if parallel
2186 vComm->AllReduce(hash, LibUtilities::ReduceSum);
2187 m_hash = hash;
2188
2189 CalculateBndSystemBandWidth();
2190 CalculateFullSystemBandWidth();
2191}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
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:77
std::set< int > m_parallelDirBndSign
Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change.
Definition: AssemblyMapCG.h:142
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Definition: AssemblyMapCG.h:106
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
Definition: AssemblyMapCG.cpp:2310
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Definition: AssemblyMapCG.h:108
void CalculateFullSystemBandWidth()
Calculate the bandwith of the full matrix system.
Definition: AssemblyMapCG.cpp:2705
std::set< ExtraDirDof > m_copyLocalDirDofs
Set indicating degrees of freedom which are Dirichlet but whose value is stored on another processor.
Definition: AssemblyMapCG.h:139
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition: AssemblyMapCG.h:112
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
Definition: AssemblyMapCG.h:130
Array< OneD, int > m_bndCondCoeffsToLocalCoeffsMap
Integer map of bnd cond coeffs to local coefficients.
Definition: AssemblyMap.h:389
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
Definition: AssemblyMap.h:387
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Definition: AssemblyMap.cpp:503
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:383
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:430
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:437
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:426
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:351
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:432
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Definition: AssemblyMap.h:381
Gs::gs_data * m_dirBndGsh
gs gather communication to impose Dirhichlet BCs.
Definition: AssemblyMap.h:422
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
Definition: AssemblyMap.h:385
Array< OneD, NekDouble > m_bndCondCoeffsToLocalCoeffsSign
Integer map of sign of bnd cond coeffs to local coefficients.
Definition: AssemblyMap.h:391
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:428
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:265
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:190
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:133
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:143
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:109
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:121
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
@ eIterativeMultiLevelStaticCond
Definition: MultiRegions.hpp:74
@ eDirectMultiLevelStaticCond
Definition: MultiRegions.hpp:71
@ eXxtMultiLevelStaticCond
Definition: MultiRegions.hpp:77
@ ePETScMultiLevelStaticCond
Definition: MultiRegions.hpp:80
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
Definition: SubStructuredGraph.h:60
std::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:51
std::vector< std::map< int, int > > DofGraph
Definition: AssemblyMapCG.h:53
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:2221
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:2269
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:50
std::vector< double > q(NPUPPER *NPUPPER)
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.hpp:644
Definition: GsLib.hpp:149
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::UnitTests::q(), Nektar::LibUtilities::ReduceSum, SetUpUniversalC0ContMap(), sign, and Vmath::Vmax().
◆ ~AssemblyMapCG()
|
override |
Destructor.
Definition at line 2196 of file AssemblyMapCG.cpp.
2197{
2200}
References Gs::Finalise(), Nektar::MultiRegions::AssemblyMap::m_bndGsh, and Nektar::MultiRegions::AssemblyMap::m_gsh.
Member Function Documentation
◆ 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 2705 of file AssemblyMapCG.cpp.
2706{
2707 int i, j;
2708 int cnt = 0;
2709 int locSize;
2710 int maxId;
2711 int minId;
2712 int bwidth = -1;
2714 {
2715 locSize =
2717 maxId = -1;
2718 minId = m_numLocalCoeffs + 1;
2719 for (j = 0; j < locSize; j++)
2720 {
2722 {
2724 {
2725 maxId = m_localToGlobalMap[cnt + j];
2726 }
2727
2729 {
2730 minId = m_localToGlobalMap[cnt + j];
2731 }
2732 }
2733 }
2734 bwidth = (bwidth > (maxId - minId)) ? bwidth : (maxId - minId);
2735
2736 cnt += locSize;
2737 }
2738
2739 m_fullSystemBandWidth = bwidth;
2740}
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
Definition: AssemblyMapCG.h:110
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()
|
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 77 of file AssemblyMapCG.cpp.
85{
86 int graphVertId = 0;
87 int vMaxVertId = -1;
88 int i, j, k, l, cnt;
89 int meshVertId, meshEdgeId, meshFaceId;
90 int meshVertId2, meshEdgeId2;
91
92 LocalRegions::ExpansionSharedPtr exp, bndExp;
95
96 m_numLocalBndCondCoeffs = 0;
97 m_systemSingular = checkIfSystemSingular;
98
99 for (i = 0; i < bndCondExp.size(); i++)
100 {
101
102 m_numLocalBndCondCoeffs += bndCondExp[i]->GetNcoeffs();
103
104 if (bndConditions[0][i]->GetBoundaryConditionType() ==
106 {
107 continue;
108 }
109
110 // Check to see if any value on boundary has Dirichlet
111 // value. note this is a vector to manage coupled
112 // solver but for scalar will just be a vector of size 11
113 cnt = 0;
114 for (k = 0; k < bndConditions.size(); ++k)
115 {
116 if (bndConditions[k][i]->GetBoundaryConditionType() ==
118 {
119 cnt++;
120 }
121 if (bndConditions[k][i]->GetBoundaryConditionType() !=
123 {
125 }
126 }
127
128 // Find the maximum boundary vertex ID on this process. This is
129 // used later to pin a vertex if the system is singular.
130 for (j = 0; j < bndCondExp[i]->GetNumElmts(); ++j)
131 {
132 bndExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::Expansion>();
133 for (k = 0; k < bndExp->GetNverts(); ++k)
134 {
135 if (vMaxVertId < bndExp->GetGeom()->GetVid(k))
136 {
137 vMaxVertId = bndExp->GetGeom()->GetVid(k);
138 }
139 }
140 }
141
142 // If all boundaries are Dirichlet fill in graph
143 if (cnt == bndConditions.size())
144 {
145 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
146 {
147 bndExp = bndCondExp[i]->GetExp(j);
148
149 for (k = 0; k < bndExp->GetNverts(); k++)
150 {
151 meshVertId = bndExp->GetGeom()->GetVid(k);
152 if (graph[0].count(meshVertId) == 0)
153 {
154 graph[0][meshVertId] = graphVertId++;
155 }
156 }
157
158 const int bndDim = bndExp->GetNumBases();
159 if (bndDim > 1)
160 {
161 for (k = 0; k < bndExp->GetNtraces(); k++)
162 {
163 meshEdgeId = bndExp->GetGeom()->GetEid(k);
164 if (graph[1].count(meshEdgeId) == 0)
165 {
166 graph[1][meshEdgeId] = graphVertId++;
167 }
168 }
169 }
170
171 // Possibility of a face in 3D or edge in 2D
172 meshFaceId = bndExp->GetGeom()->GetGlobalID();
173 if (graph[bndDim].count(meshFaceId) == 0)
174 {
175 graph[bndDim][meshFaceId] = graphVertId++;
176 }
177 m_numLocalDirBndCoeffs += bndExp->GetNcoeffs();
178 }
179 }
180 }
181
182 // Number of dirichlet edges and faces (not considering periodic
183 // BCs)
184 m_numDirEdges = graph[1].size();
185 m_numDirFaces = graph[2].size();
186
187 /*
188 * The purpose of this routine is to deal with those degrees of
189 * freedom that are Dirichlet, but do not have a local Dirichlet
190 * boundary condition expansion set.
191 *
192 * For example, in 2D, consider a triangulation of a square into two
193 * triangles. Now imagine one edge of the square is Dirichlet and
194 * the problem is run on two processors. On one processor, one
195 * triangle vertex is Dirichlet, but doesn't know this since the
196 * Dirichlet composite lives on the other processor.
197 *
198 * When the global linear system is solved therefore, there is an
199 * inconsistency that at best leads to an inaccurate answer or a
200 * divergence of the system.
201 *
202 * This routine identifies such cases for 2D, and also for 3D where
203 * e.g. edges may have the same problem (consider an extrusion of
204 * the case above, for example).
205 */
206
207 // Collate information on Dirichlet vertices from all processes
208 int n = vComm->GetSize();
210
211 if (vComm->IsSerial())
212 {
213 // for FieldConvert Comm this is true and it resets
214 // parallel processing back to serial case
215 n = 1;
216 p = 0;
217 }
218 // At this point, graph only contains information from Dirichlet
219 // boundaries. Therefore make a global list of the vert and edge
220 // information on all processors.
221 Array<OneD, int> vertcounts(n, 0);
222 Array<OneD, int> vertoffsets(n, 0);
223 Array<OneD, int> edgecounts(n, 0);
224 Array<OneD, int> edgeoffsets(n, 0);
225 vertcounts[p] = graph[0].size();
226 edgecounts[p] = graph[1].size();
227 vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
228 vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
229
230 for (i = 1; i < n; ++i)
231 {
232 vertoffsets[i] = vertoffsets[i - 1] + vertcounts[i - 1];
233 edgeoffsets[i] = edgeoffsets[i - 1] + edgecounts[i - 1];
234 }
235
238
239 Array<OneD, int> vertlist(nTotVerts, 0);
240 Array<OneD, int> edgelist(nTotEdges, 0);
241
242 // construct list of global ids of global vertices
243 i = 0;
244 for (auto &it : graph[0])
245 {
246 vertlist[vertoffsets[p] + i++] = it.first;
247 }
248
249 // construct list of global ids of global edges
250 i = 0;
251 for (auto &it : graph[1])
252 {
253 edgelist[edgeoffsets[p] + i++] = it.first;
254 }
255 vComm->AllReduce(vertlist, LibUtilities::ReduceSum);
256 vComm->AllReduce(edgelist, LibUtilities::ReduceSum);
257
258 // Now we have a list of all Dirichlet vertices and edges on all
259 // processors.
260 nExtraDirichlet = 0;
261 map<int, int> extraDirVertIds, extraDirEdgeIds;
262
263 // Ensure Dirchlet vertices are consistently recorded between
264 // processes (e.g. Dirichlet region meets Neumann region across a
265 // partition boundary requires vertex on partition to be Dirichlet).
266 //
267 // To do this we look over all elements and vertices in local
268 // partition and see if they match the values stored in the vertlist
269 // from other processors and if so record the meshVertId/meshEdgeId
270 // and the processor it comes from.
271 for (i = 0; i < n; ++i)
272 {
274 {
275 continue;
276 }
277
278 for (j = 0; j < locExpVector.size(); j++)
279 {
280 exp = locExpVector[j];
281
282 for (k = 0; k < exp->GetNverts(); k++)
283 {
284 meshVertId = exp->GetGeom()->GetVid(k);
285 if (graph[0].count(meshVertId) == 0)
286 {
287 for (l = 0; l < vertcounts[i]; ++l)
288 {
289 if (vertlist[vertoffsets[i] + l] == meshVertId)
290 {
291 extraDirVertIds[meshVertId] = i;
292 graph[0][meshVertId] = graphVertId++;
293 nExtraDirichlet++;
294 }
295 }
296 }
297 }
298
299 for (k = 0; k < exp->GetGeom()->GetNumEdges(); k++)
300 {
301 meshEdgeId = exp->GetGeom()->GetEid(k);
302 if (graph[1].count(meshEdgeId) == 0)
303 {
304 for (l = 0; l < edgecounts[i]; ++l)
305 {
306 if (edgelist[edgeoffsets[i] + l] == meshEdgeId)
307 {
308 extraDirEdgeIds[meshEdgeId] = i;
309 graph[1][meshEdgeId] = graphVertId++;
310 if (exp->GetGeom()->GetNumFaces())
311 {
312 nExtraDirichlet +=
313 exp->as<LocalRegions::Expansion3D>()
314 ->GetEdgeNcoeffs(k) -
315 2;
316 }
317 else
318 {
319 nExtraDirichlet += exp->GetTraceNcoeffs(k) - 2;
320 }
321 }
322 }
323 }
324 }
325 }
326 }
327
328 // Low Energy preconditioner needs to know how many extra Dirichlet
329 // edges are on this process so store map in array.
330 m_extraDirEdges = Array<OneD, int>(extraDirEdgeIds.size(), -1);
331 i = 0;
332 for (auto &it : extraDirEdgeIds)
333 {
334 meshEdgeId = it.first;
335 m_extraDirEdges[i++] = meshEdgeId;
336 }
337
338 // Now we have a list of all vertices and edges that are Dirichlet
339 // and not defined on the local partition as well as which processor
340 // they are stored on.
341 //
342 // Make a full list of all such entities on all processors and which
343 // processor they belong to.
344 for (i = 0; i < n; ++i)
345 {
346 vertcounts[i] = 0;
347 vertoffsets[i] = 0;
348 edgecounts[i] = 0;
349 edgeoffsets[i] = 0;
350 }
351
352 vertcounts[p] = extraDirVertIds.size();
353 edgecounts[p] = extraDirEdgeIds.size();
354 vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
355 vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
356 nTotVerts = Vmath::Vsum(n, vertcounts, 1);
357 nTotEdges = Vmath::Vsum(n, edgecounts, 1);
358
359 vertoffsets[0] = edgeoffsets[0] = 0;
360
361 for (i = 1; i < n; ++i)
362 {
363 vertoffsets[i] = vertoffsets[i - 1] + vertcounts[i - 1];
364 edgeoffsets[i] = edgeoffsets[i - 1] + edgecounts[i - 1];
365 }
366
367 Array<OneD, int> vertids(nTotVerts, 0);
368 Array<OneD, int> edgeids(nTotEdges, 0);
369 Array<OneD, int> vertprocs(nTotVerts, 0);
370 Array<OneD, int> edgeprocs(nTotEdges, 0);
371
372 i = 0;
373 for (auto &it : extraDirVertIds)
374 {
375 vertids[vertoffsets[p] + i] = it.first;
376 vertprocs[vertoffsets[p] + i] = it.second;
377 ++i;
378 }
379
380 i = 0;
381 for (auto &it : extraDirEdgeIds)
382 {
383 edgeids[edgeoffsets[p] + i] = it.first;
384 edgeprocs[edgeoffsets[p] + i] = it.second;
385 ++i;
386 }
387
388 vComm->AllReduce(vertids, LibUtilities::ReduceSum);
389 vComm->AllReduce(vertprocs, LibUtilities::ReduceSum);
390 vComm->AllReduce(edgeids, LibUtilities::ReduceSum);
391 vComm->AllReduce(edgeprocs, LibUtilities::ReduceSum);
392
393 // Set up list of vertices that need to be shared to other
394 // partitions
395 for (i = 0; i < nTotVerts; ++i)
396 {
398 {
399 extraDirVerts.insert(vertids[i]);
400 }
401 }
402
403 // Set up list of edges that need to be shared to other partitions
404 for (i = 0; i < nTotEdges; ++i)
405 {
407 {
408 extraDirEdges.insert(edgeids[i]);
409 }
410 }
411
412 // Check between processes if the whole system is singular
414 vComm->AllReduce(s, LibUtilities::ReduceMin);
416
417 // Find the minimum boundary vertex ID on each process
418 Array<OneD, int> bcminvertid(n, 0);
419 bcminvertid[p] = vMaxVertId;
420 vComm->AllReduce(bcminvertid, LibUtilities::ReduceMax);
421
422 // Find the process rank with the minimum boundary vertex ID
424
425 // If the system is singular, the process with the maximum
426 // number of BCs will set a Dirichlet vertex to make
427 // system non-singular. Note: we find the process with
428 // maximum boundary regions to ensure we do not try to set
429 // a Dirichlet vertex on a partition with no intersection
430 // with the boundary.
431 meshVertId = 0;
432
434 {
436 {
438 }
439 else if (vMaxVertId == -1)
440 {
441 // All boundaries are periodic.
442 meshVertId = locExpVector[0]->GetGeom()->GetVid(0);
443 }
444 else
445 {
446 // Set pinned vertex to that with minimum vertex ID to
447 // ensure consistency in parallel.
448 meshVertId = bcminvertid[p];
449 }
450
451 if (graph[0].count(meshVertId) == 0)
452 {
453 graph[0][meshVertId] = graphVertId++;
454 }
455 }
456
457 vComm->AllReduce(meshVertId, LibUtilities::ReduceSum);
458
459 // When running in parallel, we need to ensure that the singular
460 // mesh vertex is communicated to any periodic vertices, otherwise
461 // the system may diverge.
463 {
464 // Firstly, we check that no other processors have this
465 // vertex. If they do, then we mark the vertex as also being
466 // Dirichlet.
468 {
469 for (i = 0; i < locExpVector.size(); ++i)
470 {
471 for (j = 0; j < locExpVector[i]->GetNverts(); ++j)
472 {
473 if (locExpVector[i]->GetGeom()->GetVid(j) != meshVertId)
474 {
475 continue;
476 }
477
478 if (graph[0].count(meshVertId) == 0)
479 {
480 graph[0][meshVertId] = graphVertId++;
481 }
482 }
483 }
484 }
485
486 // In the case that meshVertId is periodic with other vertices,
487 // this process and all other processes need to make sure that
488 // the periodic vertices are also marked as Dirichlet.
489 int gId;
490
491 // At least one process (maxBCidx) will have already associated
492 // a graphVertId with meshVertId. Others won't even have any of
493 // the vertices. The logic below is designed to handle both
494 // cases.
495 if (graph[0].count(meshVertId) == 0)
496 {
497 gId = -1;
498 }
499 else
500 {
501 gId = graph[0][meshVertId];
502 }
503
504 for (auto &pIt : periodicVerts)
505 {
506 // Either the vertex is local to this processor (in which
507 // case it will be in the pIt.first position) or else
508 // meshVertId might be contained within another processor's
509 // vertex list. The if statement below covers both cases. If
510 // we find it, set as Dirichlet with the vertex id gId.
511 if (pIt.first == meshVertId)
512 {
513 gId = gId < 0 ? graphVertId++ : gId;
514 graph[0][meshVertId] = gId;
515
516 for (i = 0; i < pIt.second.size(); ++i)
517 {
518 if (pIt.second[i].isLocal)
519 {
520 graph[0][pIt.second[i].id] = graph[0][meshVertId];
521 }
522 }
523 }
524 else
525 {
526 bool found = false;
527 for (i = 0; i < pIt.second.size(); ++i)
528 {
529 if (pIt.second[i].id == meshVertId)
530 {
531 found = true;
532 break;
533 }
534 }
535
536 if (found)
537 {
538 gId = gId < 0 ? graphVertId++ : gId;
539 graph[0][pIt.first] = gId;
540
541 for (i = 0; i < pIt.second.size(); ++i)
542 {
543 if (pIt.second[i].isLocal)
544 {
545 graph[0][pIt.second[i].id] = graph[0][pIt.first];
546 }
547 }
548 }
549 }
550 }
551 }
552
553 // Add extra dirichlet boundary conditions to count.
554 m_numLocalDirBndCoeffs += nExtraDirichlet;
555 firstNonDirGraphVertId = graphVertId;
556
557 typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS>
558 BoostGraph;
559 BoostGraph boostGraphObj;
560
561 vector<map<int, int>> tempGraph(3);
562 map<int, int> vwgts_map;
563 Array<OneD, int> localVerts;
564 Array<OneD, int> localEdges;
565 Array<OneD, int> localFaces;
566
567 int tempGraphVertId = 0;
568 int localVertOffset = 0;
569 int localEdgeOffset = 0;
570 int localFaceOffset = 0;
571 int nTotalVerts = 0;
572 int nTotalEdges = 0;
573 int nTotalFaces = 0;
574 int nVerts;
575 int nEdges;
576 int nFaces;
577 int vertCnt;
578 int edgeCnt;
579 int faceCnt;
580
581 m_numNonDirVertexModes = 0;
582 m_numNonDirEdges = 0;
583 m_numNonDirFaces = 0;
584 m_numNonDirFaceModes = 0;
585 m_numNonDirFaceModes = 0;
586 m_numLocalBndCoeffs = 0;
587
588 map<int, int> EdgeSize;
589 map<int, int> FaceSize;
590
591 /// - Count verts, edges, face and add up edges and face sizes
592 for (i = 0; i < locExpVector.size(); ++i)
593 {
594 exp = locExpVector[i];
595 nEdges = exp->GetGeom()->GetNumEdges();
596 nFaces = exp->GetGeom()->GetNumFaces();
597
598 nTotalVerts += exp->GetNverts();
599 nTotalEdges += nEdges;
600 nTotalFaces += nFaces;
601
602 for (j = 0; j < nEdges; ++j)
603 {
604 meshEdgeId = exp->GetGeom()->GetEid(j);
605 int nEdgeInt;
606
607 if (nFaces)
608 {
609 nEdgeInt =
610 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
611 }
612 else
613 {
614 nEdgeInt = exp->GetTraceNcoeffs(j) - 2;
615 }
616
617 if (EdgeSize.count(meshEdgeId) > 0)
618 {
619 EdgeSize[meshEdgeId] = min(EdgeSize[meshEdgeId], nEdgeInt);
620 }
621 else
622 {
623 EdgeSize[meshEdgeId] = nEdgeInt;
624 }
625 }
626
627 faceCnt = 0;
628 for (j = 0; j < nFaces; ++j)
629 {
630 meshFaceId = exp->GetGeom()->GetFid(j);
631 if (FaceSize.count(meshFaceId) > 0)
632 {
633 FaceSize[meshFaceId] =
634 min(FaceSize[meshFaceId], exp->GetTraceIntNcoeffs(j));
635 }
636 else
637 {
638 FaceSize[meshFaceId] = exp->GetTraceIntNcoeffs(j);
639 }
640 FaceSize[meshFaceId] = exp->GetTraceIntNcoeffs(j);
641 }
642 }
643
644 /// - Periodic vertices
645 for (auto &pIt : periodicVerts)
646 {
647 meshVertId = pIt.first;
648
649 // This periodic vertex is joined to a Dirichlet condition.
650 if (graph[0].count(pIt.first) != 0)
651 {
652 for (i = 0; i < pIt.second.size(); ++i)
653 {
654 meshVertId2 = pIt.second[i].id;
655 if (graph[0].count(meshVertId2) == 0 && pIt.second[i].isLocal)
656 {
657 graph[0][meshVertId2] = graph[0][meshVertId];
658 }
659 }
660 continue;
661 }
662
663 // One of the attached vertices is Dirichlet.
664 bool isDirichlet = false;
665 for (i = 0; i < pIt.second.size(); ++i)
666 {
667 if (!pIt.second[i].isLocal)
668 {
669 continue;
670 }
671
672 meshVertId2 = pIt.second[i].id;
673 if (graph[0].count(meshVertId2) > 0)
674 {
675 isDirichlet = true;
676 break;
677 }
678 }
679
680 if (isDirichlet)
681 {
682 graph[0][meshVertId] = graph[0][pIt.second[i].id];
683
684 for (j = 0; j < pIt.second.size(); ++j)
685 {
686 meshVertId2 = pIt.second[i].id;
687 if (j == i || !pIt.second[j].isLocal ||
688 graph[0].count(meshVertId2) > 0)
689 {
690 continue;
691 }
692
693 graph[0][meshVertId2] = graph[0][pIt.second[i].id];
694 }
695
696 continue;
697 }
698
699 // Otherwise, see if a vertex ID has already been set.
700 for (i = 0; i < pIt.second.size(); ++i)
701 {
702 if (!pIt.second[i].isLocal)
703 {
704 continue;
705 }
706
707 if (tempGraph[0].count(pIt.second[i].id) > 0)
708 {
709 break;
710 }
711 }
712
713 if (i == pIt.second.size())
714 {
715 boost::add_vertex(boostGraphObj);
716 tempGraph[0][meshVertId] = tempGraphVertId++;
717 m_numNonDirVertexModes++;
718 }
719 else
720 {
721 tempGraph[0][meshVertId] = tempGraph[0][pIt.second[i].id];
722 }
723 }
724
725 // Store the temporary graph vertex id's of all element edges and
726 // vertices in these 3 arrays below
727 localVerts = Array<OneD, int>(nTotalVerts, -1);
728 localEdges = Array<OneD, int>(nTotalEdges, -1);
729 localFaces = Array<OneD, int>(nTotalFaces, -1);
730
731 // Set up vertex numbering
732 for (i = 0; i < locExpVector.size(); ++i)
733 {
734 exp = locExpVector[i];
735 vertCnt = 0;
736 nVerts = exp->GetNverts();
737 for (j = 0; j < nVerts; ++j)
738 {
739 meshVertId = exp->GetGeom()->GetVid(j);
740 if (graph[0].count(meshVertId) == 0)
741 {
742 if (tempGraph[0].count(meshVertId) == 0)
743 {
744 boost::add_vertex(boostGraphObj);
745 tempGraph[0][meshVertId] = tempGraphVertId++;
746 m_numNonDirVertexModes += 1;
747 }
748 localVerts[localVertOffset + vertCnt++] =
749 tempGraph[0][meshVertId];
750 vwgts_map[tempGraph[0][meshVertId]] = 1;
751 }
752 }
753
754 localVertOffset += nVerts;
755 }
756
757 /// - Periodic edges
758 for (auto &pIt : periodicEdges)
759 {
760 meshEdgeId = pIt.first;
761
762 // This periodic edge is joined to a Dirichlet condition.
763 if (graph[1].count(pIt.first) != 0)
764 {
765 for (i = 0; i < pIt.second.size(); ++i)
766 {
767 meshEdgeId2 = pIt.second[i].id;
768 if (graph[1].count(meshEdgeId2) == 0 && pIt.second[i].isLocal)
769 {
770 graph[1][meshEdgeId2] = graph[1][meshEdgeId];
771 }
772 }
773 continue;
774 }
775
776 // One of the attached edges is Dirichlet.
777 bool isDirichlet = false;
778 for (i = 0; i < pIt.second.size(); ++i)
779 {
780 if (!pIt.second[i].isLocal)
781 {
782 continue;
783 }
784
785 meshEdgeId2 = pIt.second[i].id;
786 if (graph[1].count(meshEdgeId2) > 0)
787 {
788 isDirichlet = true;
789 break;
790 }
791 }
792
793 if (isDirichlet)
794 {
795 graph[1][meshEdgeId] = graph[1][pIt.second[i].id];
796
797 for (j = 0; j < pIt.second.size(); ++j)
798 {
799 meshEdgeId2 = pIt.second[i].id;
800 if (j == i || !pIt.second[j].isLocal ||
801 graph[1].count(meshEdgeId2) > 0)
802 {
803 continue;
804 }
805
806 graph[1][meshEdgeId2] = graph[1][pIt.second[i].id];
807 }
808
809 continue;
810 }
811
812 // Otherwise, see if a edge ID has already been set.
813 for (i = 0; i < pIt.second.size(); ++i)
814 {
815 if (!pIt.second[i].isLocal)
816 {
817 continue;
818 }
819
820 if (tempGraph[1].count(pIt.second[i].id) > 0)
821 {
822 break;
823 }
824 }
825
826 if (i == pIt.second.size())
827 {
828 boost::add_vertex(boostGraphObj);
829 tempGraph[1][meshEdgeId] = tempGraphVertId++;
830 m_numNonDirEdgeModes += EdgeSize[meshEdgeId];
831 m_numNonDirEdges++;
832 }
833 else
834 {
835 tempGraph[1][meshEdgeId] = tempGraph[1][pIt.second[i].id];
836 }
837 }
838
839 int nEdgeIntCoeffs, nFaceIntCoeffs;
840
841 // Set up edge numbering
842 for (i = 0; i < locExpVector.size(); ++i)
843 {
844 exp = locExpVector[i];
845 edgeCnt = 0;
846 nEdges = exp->GetGeom()->GetNumEdges();
847
848 for (j = 0; j < nEdges; ++j)
849 {
850 meshEdgeId = exp->GetGeom()->GetEid(j);
851 nEdgeIntCoeffs = EdgeSize[meshEdgeId];
852 if (graph[1].count(meshEdgeId) == 0)
853 {
854 if (tempGraph[1].count(meshEdgeId) == 0)
855 {
856 boost::add_vertex(boostGraphObj);
857 tempGraph[1][meshEdgeId] = tempGraphVertId++;
858 m_numNonDirEdgeModes += nEdgeIntCoeffs;
859
860 m_numNonDirEdges++;
861 }
862 localEdges[localEdgeOffset + edgeCnt++] =
863 tempGraph[1][meshEdgeId];
864 vwgts_map[tempGraph[1][meshEdgeId]] = nEdgeIntCoeffs;
865 }
866 }
867
868 localEdgeOffset += nEdges;
869 }
870
871 /// - Periodic faces
872 for (auto &pIt : periodicFaces)
873 {
874 if (!pIt.second[0].isLocal)
875 {
876 // The face mapped to is on another process.
877 meshFaceId = pIt.first;
878 ASSERTL0(graph[2].count(meshFaceId) == 0,
879 "This periodic boundary edge has been specified before");
880 boost::add_vertex(boostGraphObj);
881 tempGraph[2][meshFaceId] = tempGraphVertId++;
882 nFaceIntCoeffs = FaceSize[meshFaceId];
883 m_numNonDirFaceModes += nFaceIntCoeffs;
884 m_numNonDirFaces++;
885 }
886 else if (pIt.first < pIt.second[0].id)
887 {
888 ASSERTL0(graph[2].count(pIt.first) == 0,
889 "This periodic boundary face has been specified before");
890 ASSERTL0(graph[2].count(pIt.second[0].id) == 0,
891 "This periodic boundary face has been specified before");
892
893 boost::add_vertex(boostGraphObj);
894 tempGraph[2][pIt.first] = tempGraphVertId;
895 tempGraph[2][pIt.second[0].id] = tempGraphVertId++;
896 nFaceIntCoeffs = FaceSize[pIt.first];
897 m_numNonDirFaceModes += nFaceIntCoeffs;
898 m_numNonDirFaces++;
899 }
900 }
901
902 // setup face numbering
903 for (i = 0; i < locExpVector.size(); ++i)
904 {
905 exp = locExpVector[i];
906 nFaces = exp->GetGeom()->GetNumFaces();
907 faceCnt = 0;
908 for (j = 0; j < nFaces; ++j)
909 {
910 nFaceIntCoeffs = exp->GetTraceIntNcoeffs(j);
911 meshFaceId = exp->GetGeom()->GetFid(j);
912 if (graph[2].count(meshFaceId) == 0)
913 {
914 if (tempGraph[2].count(meshFaceId) == 0)
915 {
916 boost::add_vertex(boostGraphObj);
917 tempGraph[2][meshFaceId] = tempGraphVertId++;
918 m_numNonDirFaceModes += nFaceIntCoeffs;
919
920 m_numNonDirFaces++;
921 }
922 localFaces[localFaceOffset + faceCnt++] =
923 tempGraph[2][meshFaceId];
924 vwgts_map[tempGraph[2][meshFaceId]] = nFaceIntCoeffs;
925 }
926 }
927 m_numLocalBndCoeffs += exp->NumBndryCoeffs();
928
929 localFaceOffset += nFaces;
930 }
931
932 localVertOffset = 0;
933 localEdgeOffset = 0;
934 localFaceOffset = 0;
935 for (i = 0; i < locExpVector.size(); ++i)
936 {
937 exp = locExpVector[i];
938 nVerts = exp->GetNverts();
939 nEdges = exp->GetGeom()->GetNumEdges();
940 nFaces = exp->GetGeom()->GetNumFaces();
941
942 // Now loop over all local faces, edges and vertices of this
943 // element and define that all other faces, edges and verices of
944 // this element are adjacent to them.
945
946 // Vertices
947 for (j = 0; j < nVerts; j++)
948 {
949 if (localVerts[j + localVertOffset] == -1)
950 {
951 break;
952 }
953 // associate to other vertices
954 for (k = 0; k < nVerts; k++)
955 {
956 if (localVerts[k + localVertOffset] == -1)
957 {
958 break;
959 }
960 if (k != j)
961 {
962 boost::add_edge((size_t)localVerts[j + localVertOffset],
963 (size_t)localVerts[k + localVertOffset],
964 boostGraphObj);
965 }
966 }
967 // associate to other edges
968 for (k = 0; k < nEdges; k++)
969 {
970 if (localEdges[k + localEdgeOffset] == -1)
971 {
972 break;
973 }
974 boost::add_edge((size_t)localVerts[j + localVertOffset],
975 (size_t)localEdges[k + localEdgeOffset],
976 boostGraphObj);
977 }
978 // associate to other faces
979 for (k = 0; k < nFaces; k++)
980 {
981 if (localFaces[k + localFaceOffset] == -1)
982 {
983 break;
984 }
985 boost::add_edge((size_t)localVerts[j + localVertOffset],
986 (size_t)localFaces[k + localFaceOffset],
987 boostGraphObj);
988 }
989 }
990
991 // Edges
992 for (j = 0; j < nEdges; j++)
993 {
994 if (localEdges[j + localEdgeOffset] == -1)
995 {
996 break;
997 }
998 // Associate to other edges
999 for (k = 0; k < nEdges; k++)
1000 {
1001 if (localEdges[k + localEdgeOffset] == -1)
1002 {
1003 break;
1004 }
1005 if (k != j)
1006 {
1007 boost::add_edge((size_t)localEdges[j + localEdgeOffset],
1008 (size_t)localEdges[k + localEdgeOffset],
1009 boostGraphObj);
1010 }
1011 }
1012 // Associate to vertices
1013 for (k = 0; k < nVerts; k++)
1014 {
1015 if (localVerts[k + localVertOffset] == -1)
1016 {
1017 break;
1018 }
1019 boost::add_edge((size_t)localEdges[j + localEdgeOffset],
1020 (size_t)localVerts[k + localVertOffset],
1021 boostGraphObj);
1022 }
1023 // Associate to faces
1024 for (k = 0; k < nFaces; k++)
1025 {
1026 if (localFaces[k + localFaceOffset] == -1)
1027 {
1028 break;
1029 }
1030 boost::add_edge((size_t)localEdges[j + localEdgeOffset],
1031 (size_t)localFaces[k + localFaceOffset],
1032 boostGraphObj);
1033 }
1034 }
1035
1036 // Faces
1037 for (j = 0; j < nFaces; j++)
1038 {
1039 if (localFaces[j + localFaceOffset] == -1)
1040 {
1041 break;
1042 }
1043 // Associate to other faces
1044 for (k = 0; k < nFaces; k++)
1045 {
1046 if (localFaces[k + localFaceOffset] == -1)
1047 {
1048 break;
1049 }
1050 if (k != j)
1051 {
1052 boost::add_edge((size_t)localFaces[j + localFaceOffset],
1053 (size_t)localFaces[k + localFaceOffset],
1054 boostGraphObj);
1055 }
1056 }
1057 // Associate to vertices
1058 for (k = 0; k < nVerts; k++)
1059 {
1060 if (localVerts[k + localVertOffset] == -1)
1061 {
1062 break;
1063 }
1064 boost::add_edge((size_t)localFaces[j + localFaceOffset],
1065 (size_t)localVerts[k + localVertOffset],
1066 boostGraphObj);
1067 }
1068 // Associate to edges
1069 for (k = 0; k < nEdges; k++)
1070 {
1071 if (localEdges[k + localEdgeOffset] == -1)
1072 {
1073 break;
1074 }
1075 boost::add_edge((size_t)localFaces[j + localFaceOffset],
1076 (size_t)localEdges[k + localEdgeOffset],
1077 boostGraphObj);
1078 }
1079 }
1080
1081 localVertOffset += nVerts;
1082 localEdgeOffset += nEdges;
1083 localFaceOffset += nFaces;
1084 }
1085
1086 // Container to store vertices of the graph which correspond to
1087 // degrees of freedom along the boundary and periodic BCs.
1088 set<int> partVerts;
1089
1092 {
1093 vector<long> procVerts, procEdges, procFaces;
1094 set<int> foundVerts, foundEdges, foundFaces;
1095
1096 // Loop over element and construct the procVerts and procEdges
1097 // vectors, which store the geometry IDs of mesh vertices and
1098 // edges respectively which are local to this process.
1099 for (i = cnt = 0; i < locExpVector.size(); ++i)
1100 {
1101 int elmtid = i;
1102 exp = locExpVector[elmtid];
1103 for (j = 0; j < exp->GetNverts(); ++j)
1104 {
1105 int vid = exp->GetGeom()->GetVid(j) + 1;
1106 if (foundVerts.count(vid) == 0)
1107 {
1108 procVerts.push_back(vid);
1109 foundVerts.insert(vid);
1110 }
1111 }
1112
1113 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
1114 {
1115 int eid = exp->GetGeom()->GetEid(j) + 1;
1116
1117 if (foundEdges.count(eid) == 0)
1118 {
1119 procEdges.push_back(eid);
1120 foundEdges.insert(eid);
1121 }
1122 }
1123
1124 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
1125 {
1126 int fid = exp->GetGeom()->GetFid(j) + 1;
1127
1128 if (foundFaces.count(fid) == 0)
1129 {
1130 procFaces.push_back(fid);
1131 foundFaces.insert(fid);
1132 }
1133 }
1134 }
1135
1136 int unique_verts = foundVerts.size();
1137 int unique_edges = foundEdges.size();
1138 int unique_faces = foundFaces.size();
1139
1141
1142 // Now construct temporary GS objects. These will be used to
1143 // populate the arrays tmp3 and tmp4 with the multiplicity of
1144 // the vertices and edges respectively to identify those
1145 // vertices and edges which are located on partition boundary.
1146 Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1148 Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1149 Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1150 Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1152 Gs::Finalise(tmp1);
1153
1154 if (unique_edges > 0)
1155 {
1156 Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1159 Gs::Finalise(tmp2);
1160 }
1161
1162 if (unique_faces > 0)
1163 {
1164 Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1167 Gs::Finalise(tmp3);
1168 }
1169
1170 // Finally, fill the partVerts set with all non-Dirichlet
1171 // vertices which lie on a partition boundary.
1172 for (i = 0; i < unique_verts; ++i)
1173 {
1174 if (tmp4[i] > 1.0)
1175 {
1176 if (graph[0].count(procVerts[i] - 1) == 0)
1177 {
1178 partVerts.insert(tempGraph[0][procVerts[i] - 1]);
1179 }
1180 }
1181 }
1182
1183 for (i = 0; i < unique_edges; ++i)
1184 {
1185 if (tmp5[i] > 1.0)
1186 {
1187 if (graph[1].count(procEdges[i] - 1) == 0)
1188 {
1189 partVerts.insert(tempGraph[1][procEdges[i] - 1]);
1190 }
1191 }
1192 }
1193
1194 for (i = 0; i < unique_faces; ++i)
1195 {
1196 if (tmp6[i] > 1.0)
1197 {
1198 if (graph[2].count(procFaces[i] - 1) == 0)
1199 {
1200 partVerts.insert(tempGraph[2][procFaces[i] - 1]);
1201 }
1202 }
1203 }
1204
1205 // Now fill with all vertices on periodic BCs
1206 for (auto &pIt : periodicVerts)
1207 {
1208 if (graph[0].count(pIt.first) == 0)
1209 {
1210 partVerts.insert(tempGraph[0][pIt.first]);
1211 }
1212 }
1213 for (auto &pIt : periodicEdges)
1214 {
1215 if (graph[1].count(pIt.first) == 0)
1216 {
1217 partVerts.insert(tempGraph[1][pIt.first]);
1218 }
1219 }
1220 for (auto &pIt : periodicFaces)
1221 {
1222 if (graph[2].count(pIt.first) == 0)
1223 {
1224 partVerts.insert(tempGraph[2][pIt.first]);
1225 }
1226 }
1227 }
1228
1229 int nGraphVerts = tempGraphVertId;
1230 Array<OneD, int> perm(nGraphVerts);
1231 Array<OneD, int> iperm(nGraphVerts);
1232 Array<OneD, int> vwgts(nGraphVerts);
1234 for (i = 0; i < nGraphVerts; ++i)
1235 {
1236 vwgts[i] = vwgts_map[i];
1237 }
1238
1239 if (nGraphVerts)
1240 {
1242 {
1250 {
1251 NoReordering(boostGraphObj, perm, iperm);
1252 break;
1253 }
1254
1256 {
1257 CuthillMckeeReordering(boostGraphObj, perm, iperm);
1258 break;
1259 }
1260
1265 {
1266 MultiLevelBisectionReordering(boostGraphObj, perm, iperm,
1267 bottomUpGraph, partVerts,
1268 mdswitch);
1269 break;
1270 }
1271 default:
1272 {
1274 "Unrecognised solution type: " +
1276 }
1277 }
1278 }
1279
1280 // For parallel multi-level static condensation determine the lowest
1281 // static condensation level amongst processors.
1286 bottomUpGraph)
1287 {
1288 m_lowestStaticCondLevel = bottomUpGraph->GetNlevels() - 1;
1290 }
1291 else
1292 {
1293 m_lowestStaticCondLevel = 0;
1294 }
1295
1296 /**
1297 * STEP 4: Fill the #graph[0] and
1298 * #graph[1] with the optimal ordering from boost.
1299 */
1300 for (auto &mapIt : tempGraph[0])
1301 {
1302 graph[0][mapIt.first] = iperm[mapIt.second] + graphVertId;
1303 }
1304 for (auto &mapIt : tempGraph[1])
1305 {
1306 graph[1][mapIt.first] = iperm[mapIt.second] + graphVertId;
1307 }
1308 for (auto &mapIt : tempGraph[2])
1309 {
1310 graph[2][mapIt.first] = iperm[mapIt.second] + graphVertId;
1311 }
1312
1313 return nGraphVerts;
1314}
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:349
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
Definition: SubStructuredGraph.cpp:1012
void CuthillMckeeReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
Definition: SubStructuredGraph.cpp:682
void MultiLevelBisectionReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm, BottomUpSubStructuredGraphSharedPtr &substructgraph, std::set< int > partVerts, int mdswitch)
Definition: SubStructuredGraph.cpp:703
const char *const GlobalSysSolnTypeMap[]
Definition: MultiRegions.hpp:84
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.hpp:623
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()
|
inline |
Definition at line 94 of file AssemblyMapCG.h.
References m_copyLocalDirDofs.
◆ GetParallelDirBndSign()
|
inline |
Definition at line 99 of file AssemblyMapCG.h.
References m_parallelDirBndSign.
◆ SetUpUniversalC0ContMap()
|
protected |
Sets up the global to universal mapping of degrees of freedom across processors.
Definition at line 2310 of file AssemblyMapCG.cpp.
2314{
2315 LocalRegions::ExpansionSharedPtr exp;
2316 int nVert = 0;
2317 int nEdge = 0;
2318 int nFace = 0;
2319 int maxEdgeDof = 0;
2320 int maxFaceDof = 0;
2321 int maxIntDof = 0;
2322 int dof = 0;
2323 int cnt;
2324 int i, j, k, l;
2325 int meshVertId;
2326 int meshEdgeId;
2327 int meshFaceId;
2328 int elementId;
2329 int vGlobalId;
2330 int maxBndGlobalId = 0;
2331 StdRegions::Orientation edgeOrient;
2332 StdRegions::Orientation faceOrient;
2333 Array<OneD, unsigned int> edgeInteriorMap;
2334 Array<OneD, int> edgeInteriorSign;
2335 Array<OneD, unsigned int> faceInteriorMap;
2336 Array<OneD, int> faceInteriorSign;
2337 Array<OneD, unsigned int> interiorMap;
2338
2341 const bool verbose = locExp.GetSession()->DefinesCmdLineArgument("verbose");
2342
2344 m_globalToUniversalMapUnique =
2346 m_globalToUniversalBndMap =
2350
2351 // Loop over all the elements in the domain to gather mesh data
2352 for (i = 0; i < locExpVector.size(); ++i)
2353 {
2354 exp = locExpVector[i];
2355
2356 int nv = exp->GetNverts();
2357 int ne = exp->GetGeom()->GetNumEdges();
2358 int nf = exp->GetGeom()->GetNumFaces();
2359
2360 nVert += nv;
2361 nEdge += ne;
2362 nFace += nf;
2363
2364 // Loop over all edges (and vertices) of element i
2365 for (j = 0; j < ne; ++j)
2366 {
2367 if (nf)
2368 {
2369 dof =
2370 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
2371 }
2372 else
2373 {
2374 dof = exp->GetTraceNcoeffs(j) - 2;
2375 }
2376
2377 maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2378 }
2379 for (j = 0; j < nf; ++j)
2380 {
2381 dof = exp->GetTraceIntNcoeffs(j);
2382 maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2383 }
2384 exp->GetInteriorMap(interiorMap);
2385 dof = interiorMap.size();
2386 maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2387 }
2388
2389 // Tell other processes about how many dof we have
2390 vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2391 vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2392 vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2393 vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2394 vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2395 vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2396
2397 // Assemble global to universal mapping for this process
2398 for (i = 0; i < locExpVector.size(); ++i)
2399 {
2400 exp = locExpVector[i];
2401 cnt = locExp.GetCoeff_Offset(i);
2402
2403 int nf = exp->GetGeom()->GetNumFaces();
2404
2405 // Loop over all vertices of element i
2406 for (j = 0; j < exp->GetNverts(); ++j)
2407 {
2408 meshVertId = exp->GetGeom()->GetVid(j);
2409 vGlobalId = m_localToGlobalMap[cnt + exp->GetVertexMap(j)];
2410
2411 auto pIt = perVerts.find(meshVertId);
2412 if (pIt != perVerts.end())
2413 {
2414 for (k = 0; k < pIt->second.size(); ++k)
2415 {
2416 meshVertId = min(meshVertId, pIt->second[k].id);
2417 }
2418 }
2419
2420 m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2421 m_globalToUniversalBndMap[vGlobalId] =
2422 m_globalToUniversalMap[vGlobalId];
2423 maxBndGlobalId =
2424 (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2425 }
2426
2427 // Loop over all edges of element i
2428 for (j = 0; j < exp->GetGeom()->GetNumEdges(); ++j)
2429 {
2430 meshEdgeId = exp->GetGeom()->GetEid(j);
2431 auto pIt = perEdges.find(meshEdgeId);
2432 edgeOrient = exp->GetGeom()->GetEorient(j);
2433
2434 if (pIt != perEdges.end())
2435 {
2436 pair<int, StdRegions::Orientation> idOrient =
2437 DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
2438 pIt->second);
2439 meshEdgeId = idOrient.first;
2440 edgeOrient = idOrient.second;
2441 }
2442
2443 if (nf) // 3D version
2444 {
2445 exp->as<LocalRegions::Expansion3D>()
2446 ->GetEdgeInteriorToElementMap(j, edgeInteriorMap,
2447 edgeInteriorSign, edgeOrient);
2448 dof =
2449 exp->as<LocalRegions::Expansion3D>()->GetEdgeNcoeffs(j) - 2;
2450 }
2451 else // 2D version
2452 {
2453 exp->GetTraceInteriorToElementMap(j, edgeInteriorMap,
2454 edgeInteriorSign, edgeOrient);
2455 dof = exp->GetTraceNcoeffs(j) - 2;
2456 }
2457
2458 // Set the global DOF's for the interior modes of edge j
2459 // for varP, ignore modes with sign == 0
2460 for (k = 0, l = 0; k < dof; ++k)
2461 {
2463 {
2465 {
2466 continue;
2467 }
2468 }
2469 vGlobalId = m_localToGlobalMap[cnt + edgeInteriorMap[k]];
2470 m_globalToUniversalMap[vGlobalId] =
2471 nVert + meshEdgeId * maxEdgeDof + l + 1;
2472 m_globalToUniversalBndMap[vGlobalId] =
2473 m_globalToUniversalMap[vGlobalId];
2474 maxBndGlobalId =
2475 (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2476 l++;
2477 }
2478 }
2479
2480 // Loop over all faces of element i
2481 for (j = 0; j < exp->GetGeom()->GetNumFaces(); ++j)
2482 {
2483 faceOrient = exp->GetGeom()->GetForient(j);
2484
2485 meshFaceId = exp->GetGeom()->GetFid(j);
2486
2487 auto pIt = perFaces.find(meshFaceId);
2488 if (pIt != perFaces.end())
2489 {
2490 if (meshFaceId == min(meshFaceId, pIt->second[0].id))
2491 {
2492 faceOrient = DeterminePeriodicFaceOrient(
2493 faceOrient, pIt->second[0].orient);
2494 }
2495 meshFaceId = min(meshFaceId, pIt->second[0].id);
2496 }
2497
2498 exp->GetTraceInteriorToElementMap(j, faceInteriorMap,
2499 faceInteriorSign, faceOrient);
2500 dof = exp->GetTraceIntNcoeffs(j);
2501
2502 for (k = 0, l = 0; k < dof; ++k)
2503 {
2505 {
2507 {
2508 continue;
2509 }
2510 }
2511 vGlobalId = m_localToGlobalMap[cnt + faceInteriorMap[k]];
2512 m_globalToUniversalMap[vGlobalId] = nVert + nEdge * maxEdgeDof +
2513 meshFaceId * maxFaceDof +
2514 l + 1;
2515 m_globalToUniversalBndMap[vGlobalId] =
2516 m_globalToUniversalMap[vGlobalId];
2517
2518 maxBndGlobalId =
2519 (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2520 l++;
2521 }
2522 }
2523
2524 // Add interior DOFs to complete universal numbering
2525 exp->GetInteriorMap(interiorMap);
2526 dof = interiorMap.size();
2527 elementId = (exp->GetGeom())->GetGlobalID();
2528 for (k = 0; k < dof; ++k)
2529 {
2530 vGlobalId = m_localToGlobalMap[cnt + interiorMap[k]];
2531 m_globalToUniversalMap[vGlobalId] = nVert + nEdge * maxEdgeDof +
2532 nFace * maxFaceDof +
2533 elementId * maxIntDof + k + 1;
2534 }
2535 }
2536
2537 // Set up the GSLib universal assemble mapping
2538 // Internal DOF do not participate in any data
2539 // exchange, so we keep these set to the special GSLib id=0 so
2540 // they are ignored.
2545 {
2546 tmp[i] = m_globalToUniversalMap[i];
2547 }
2548
2551 Gs::Unique(tmp, vCommRow);
2553 {
2554 m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2555 }
2557 {
2558 m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2559 }
2560}
Definition: SharedArray.hpp:51
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMapCG.h:114
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:397
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:399
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:225
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]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2860 of file AssemblyMapCG.cpp.
2862{
2863 Array<OneD, const NekDouble> local;
2865 {
2867 }
2868 else
2869 {
2871 }
2872
2874
2876 {
2878 m_localToGlobalMap.get(), global.get());
2879 }
2880 else
2881 {
2883 global.get());
2884 }
2885 UniversalAssemble(global);
2886}
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Definition: AssemblyMap.cpp:816
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.hpp:577
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]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2888 of file AssemblyMapCG.cpp.
2890{
2892}
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Definition: AssemblyMap.cpp:804
References Nektar::MultiRegions::AssemblyMap::Assemble(), Nektar::NekVector< DataType >::GetPtr(), and CG_Iterations::loc.
◆ v_GetExtraDirEdges()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2953 of file AssemblyMapCG.cpp.
References m_extraDirEdges.
◆ v_GetFullSystemBandWidth()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2913 of file AssemblyMapCG.cpp.
References m_fullSystemBandWidth.
◆ v_GetGlobalToUniversalMap() [1/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2762 of file AssemblyMapCG.cpp.
References m_globalToUniversalMap.
◆ v_GetGlobalToUniversalMap() [2/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2747 of file AssemblyMapCG.cpp.
References m_globalToUniversalMap.
◆ v_GetGlobalToUniversalMapUnique() [1/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2767 of file AssemblyMapCG.cpp.
References m_globalToUniversalMapUnique.
◆ v_GetGlobalToUniversalMapUnique() [2/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2752 of file AssemblyMapCG.cpp.
References m_globalToUniversalMapUnique.
◆ v_GetLocalToGlobalMap() [1/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2757 of file AssemblyMapCG.cpp.
References m_localToGlobalMap.
◆ v_GetLocalToGlobalMap() [2/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2742 of file AssemblyMapCG.cpp.
References m_localToGlobalMap.
◆ v_GetLocalToGlobalSign() [1/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2785 of file AssemblyMapCG.cpp.
References m_localToGlobalSign.
◆ v_GetLocalToGlobalSign() [2/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2773 of file AssemblyMapCG.cpp.
2774{
2776 {
2778 }
2779 else
2780 {
2781 return 1.0;
2782 }
2783}
References m_localToGlobalSign, and Nektar::MultiRegions::AssemblyMap::m_signChange.
◆ v_GetNumDirEdges()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2933 of file AssemblyMapCG.cpp.
References m_numDirEdges.
◆ v_GetNumDirFaces()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2938 of file AssemblyMapCG.cpp.
References m_numDirFaces.
◆ v_GetNumNonDirEdgeModes()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2923 of file AssemblyMapCG.cpp.
References m_numNonDirEdgeModes.
◆ v_GetNumNonDirEdges()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2943 of file AssemblyMapCG.cpp.
References m_numNonDirEdges.
◆ v_GetNumNonDirFaceModes()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2928 of file AssemblyMapCG.cpp.
References m_numNonDirFaceModes.
◆ v_GetNumNonDirFaces()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2948 of file AssemblyMapCG.cpp.
References m_numNonDirFaces.
◆ v_GetNumNonDirVertexModes()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2918 of file AssemblyMapCG.cpp.
References m_numNonDirVertexModes.
◆ v_GlobalToLocal() [1/2]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2829 of file AssemblyMapCG.cpp.
2831{
2832 Array<OneD, const NekDouble> glo;
2834 {
2835 glo = Array<OneD, NekDouble>(m_numGlobalCoeffs, global.data());
2836 }
2837 else
2838 {
2839 glo = global; // create reference
2840 }
2841
2843 {
2846 }
2847 else
2848 {
2850 loc.get());
2851 }
2852}
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]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2854 of file AssemblyMapCG.cpp.
2856{
2858}
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
Definition: AssemblyMap.cpp:792
References Nektar::NekVector< DataType >::GetPtr(), Nektar::MultiRegions::AssemblyMap::GlobalToLocal(), and CG_Iterations::loc.
◆ v_LinearSpaceMap()
|
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 2570 of file AssemblyMapCG.cpp.
2572{
2573 AssemblyMapCGSharedPtr returnval;
2574
2575 int i, j;
2576 int nverts = 0;
2577 const std::shared_ptr<LocalRegions::ExpansionVector> exp = locexp.GetExp();
2578 int nelmts = exp->size();
2579 const bool verbose = locexp.GetSession()->DefinesCmdLineArgument("verbose");
2580
2581 // Get Default Map and turn off any searched values.
2582 returnval = MemoryManager<AssemblyMapCG>::AllocateSharedPtr(
2583 m_session, locexp.GetComm());
2584 returnval->m_solnType = solnType;
2585 returnval->m_preconType = "Null";
2586 returnval->m_maxStaticCondLevel = 0;
2587 returnval->m_signChange = false;
2588 returnval->m_comm = m_comm;
2589
2590 // Count the number of vertices
2591 for (i = 0; i < nelmts; ++i)
2592 {
2593 nverts += (*exp)[i]->GetNverts();
2594 }
2595
2596 returnval->m_numLocalCoeffs = nverts;
2597 returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2598
2599 // Store original global ids in this map
2600 returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2601
2602 int cnt = 0;
2603 int cnt1 = 0;
2604 Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2605
2606 // Set up local to global map;
2607 for (i = 0; i < nelmts; ++i)
2608 {
2609 for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2610 {
2611 returnval->m_localToGlobalMap[cnt] =
2612 returnval->m_localToGlobalBndMap[cnt] =
2614 GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2615
2616 // Set up numLocalDirBndCoeffs
2618 {
2619 returnval->m_numLocalDirBndCoeffs++;
2620 }
2621 cnt++;
2622 }
2623 cnt1 += (*exp)[i]->GetNcoeffs();
2624 }
2625
2626 cnt = 0;
2627 // Reset global numbering and count number of dofs
2629 {
2630 if (GlobCoeffs[i] != -1)
2631 {
2632 GlobCoeffs[i] = cnt++;
2633 }
2634 }
2635
2636 // Set up number of globalCoeffs;
2637 returnval->m_numGlobalCoeffs = cnt;
2638
2639 // Set up number of global Dirichlet boundary coefficients
2641 {
2642 if (GlobCoeffs[i] != -1)
2643 {
2644 returnval->m_numGlobalDirBndCoeffs++;
2645 }
2646 }
2647
2648 // Set up global to universal map
2650 {
2651 LibUtilities::CommSharedPtr vCommRow =
2652 m_session->GetComm()->GetRowComm();
2653 int nglocoeffs = returnval->m_numGlobalCoeffs;
2654 returnval->m_globalToUniversalMap = Array<OneD, int>(nglocoeffs);
2655 returnval->m_globalToUniversalMapUnique = Array<OneD, int>(nglocoeffs);
2656
2657 // Reset local to global map and setup universal map
2658 for (i = 0; i < nverts; ++i)
2659 {
2660 cnt = returnval->m_localToGlobalMap[i];
2661 returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2662
2663 returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2664 m_globalToUniversalMap[cnt];
2665 }
2666
2667 Nektar::Array<OneD, long> tmp(nglocoeffs);
2668 Vmath::Zero(nglocoeffs, tmp, 1);
2669 for (unsigned int i = 0; i < nglocoeffs; ++i)
2670 {
2671 tmp[i] = returnval->m_globalToUniversalMap[i];
2672 }
2673 returnval->m_gsh = Gs::Init(tmp, vCommRow, verbose);
2674 Gs::Unique(tmp, vCommRow);
2675 for (unsigned int i = 0; i < nglocoeffs; ++i)
2676 {
2677 returnval->m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2678 }
2679 }
2680 else // not sure this option is ever needed.
2681 {
2682 for (i = 0; i < nverts; ++i)
2683 {
2684 cnt = returnval->m_localToGlobalMap[i];
2685 returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2686 }
2687 }
2688
2689 return returnval;
2690}
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:50
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), 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]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2790 of file AssemblyMapCG.cpp.
2793{
2794 Array<OneD, const NekDouble> local;
2796 {
2798 }
2799 else
2800 {
2802 }
2803
2805 {
2807 m_localToGlobalMap.get(), global.get());
2808 }
2809 else
2810 {
2812 global.get());
2813 }
2814
2815 // ensure all values are unique by calling a max
2816 if (useComm)
2817 {
2819 }
2820}
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.hpp:539
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]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2822 of file AssemblyMapCG.cpp.
2825{
2827}
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Definition: AssemblyMap.cpp:778
References Nektar::NekVector< DataType >::GetPtr(), CG_Iterations::loc, and Nektar::MultiRegions::AssemblyMap::LocalToGlobal().
◆ v_UniversalAssemble() [1/3]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2894 of file AssemblyMapCG.cpp.
References Gs::Gather(), Gs::gs_add, and Nektar::MultiRegions::AssemblyMap::m_gsh.
◆ v_UniversalAssemble() [2/3]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2904 of file AssemblyMapCG.cpp.
2906{
2907 Array<OneD, NekDouble> tmp(offset);
2908 Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2909 UniversalAssemble(pGlobal);
2910 Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2911}
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
References Nektar::MultiRegions::AssemblyMap::UniversalAssemble(), and Vmath::Vcopy().
◆ v_UniversalAssemble() [3/3]
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::AssemblyMap.
Definition at line 2899 of file AssemblyMapCG.cpp.
References Nektar::NekVector< DataType >::GetPtr(), and Nektar::MultiRegions::AssemblyMap::UniversalAssemble().
Member Data Documentation
◆ m_copyLocalDirDofs
|
protected |
Set indicating degrees of freedom which are Dirichlet but whose value is stored on another processor.
Definition at line 139 of file AssemblyMapCG.h.
Referenced by AssemblyMapCG(), and GetCopyLocalDirDofs().
◆ m_extraDirEdges
Extra dirichlet edges in parallel.
Definition at line 132 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetExtraDirEdges().
◆ m_fullSystemBandWidth
|
protected |
Bandwith of the full matrix system (no static condensation).
Definition at line 110 of file AssemblyMapCG.h.
Referenced by CalculateFullSystemBandWidth(), and v_GetFullSystemBandWidth().
◆ m_globalToUniversalMap
Integer map of process coeffs to universal space.
Definition at line 112 of file AssemblyMapCG.h.
Referenced by AssemblyMapCG(), Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), SetUpUniversalC0ContMap(), v_GetGlobalToUniversalMap(), and v_LinearSpaceMap().
◆ m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition at line 114 of file AssemblyMapCG.h.
Referenced by Nektar::CoupledAssemblyMap::CoupledAssemblyMap(), SetUpUniversalC0ContMap(), and v_GetGlobalToUniversalMapUnique().
◆ m_localToGlobalMap
Integer map of local coeffs to global space.
Definition at line 106 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
Integer sign of local coeffs to global space.
Definition at line 108 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
|
protected |
Maximum static condensation level.
Definition at line 136 of file AssemblyMapCG.h.
Referenced by AssemblyMapCG().
◆ m_numDirEdges
|
protected |
Number of Dirichlet edges.
Definition at line 122 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumDirEdges().
◆ m_numDirFaces
|
protected |
Number of Dirichlet faces.
Definition at line 124 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumDirFaces().
◆ m_numLocalBndCondCoeffs
|
protected |
Number of local boundary condition coefficients.
Definition at line 130 of file AssemblyMapCG.h.
Referenced by AssemblyMapCG(), and CreateGraph().
◆ m_numLocDirBndCondDofs
|
protected |
Number of local boundary condition degrees of freedom.
Definition at line 134 of file AssemblyMapCG.h.
◆ m_numNonDirEdgeModes
|
protected |
Number of non Dirichlet edge modes.
Definition at line 118 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumNonDirEdgeModes().
◆ m_numNonDirEdges
|
protected |
Number of Dirichlet edges.
Definition at line 126 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumNonDirEdges().
◆ m_numNonDirFaceModes
|
protected |
Number of non Dirichlet face modes.
Definition at line 120 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumNonDirFaceModes().
◆ m_numNonDirFaces
|
protected |
Number of Dirichlet faces.
Definition at line 128 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumNonDirFaces().
◆ m_numNonDirVertexModes
|
protected |
Number of non Dirichlet vertex modes.
Definition at line 116 of file AssemblyMapCG.h.
Referenced by CreateGraph(), and v_GetNumNonDirVertexModes().
◆ m_parallelDirBndSign
|
protected |
Set indicating the local coeffs just touching parallel dirichlet boundary that have a sign change.
Definition at line 142 of file AssemblyMapCG.h.
Referenced by AssemblyMapCG(), and GetParallelDirBndSign().
