#include <CoupledLocalToGlobalC0ContMap.h>

Inheritance diagram for Nektar::CoupledLocalToGlobalC0ContMap:

Collaboration diagram for Nektar::CoupledLocalToGlobalC0ContMap:

Public Member Functions
	CoupledLocalToGlobalC0ContMap (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const MultiRegions::ExpListSharedPtr &pressure, const int nz_loc, const bool CheeckForSingularSys=true)

void	FindEdgeIdToAddMeanPressure (vector< map< int, int > > &ReorderedGraphVertId, int &nel, const LocalRegions::ExpansionVector &locExpVector, int &edgeId, int &vertId, int &firstNonDirGraphVertId, map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)

Public Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
	AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
	Default constructor. More...

	AssemblyMapCG (const LibUtilities::SessionReaderSharedPtr &pSession, const int numLocalCoeffs, const ExpList &locExp, const BndCondExp &bndCondExp=NullExpListSharedPtrArray, const BndCond &bndConditions=SpatialDomains::NullBoundaryConditionShPtrArray, const bool checkIfSingular=false, const std::string variable="defaultVar", const PeriodicMap &periodicVerts=NullPeriodicMap, const PeriodicMap &periodicEdges=NullPeriodicMap, const PeriodicMap &periodicFaces=NullPeriodicMap)
	General constructor for expansions of all dimensions without boundary conditions. More...

virtual	~AssemblyMapCG ()
	Destructor. More...

map< int, vector< ExtraDirDof > > &	GetExtraDirDofs ()

Public Member Functions inherited from Nektar::MultiRegions::AssemblyMap
	AssemblyMap ()
	Default constructor. More...

	AssemblyMap (const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
	Constructor with a communicator. More...

	AssemblyMap (AssemblyMap *oldLevelMap, const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
	Constructor for next level in multi-level static condensation. More...

virtual	~AssemblyMap ()
	Destructor. More...

LibUtilities::CommSharedPtr	GetComm ()
	Retrieves the communicator. More...

size_t	GetHash () const
	Retrieves the hash of this map. More...

int	GetLocalToGlobalMap (const int i) const

int	GetGlobalToUniversalMap (const int i) const

int	GetGlobalToUniversalMapUnique (const int i) const

const Array< OneD, const int > &	GetLocalToGlobalMap ()

const Array< OneD, const int > &	GetGlobalToUniversalMap ()

const Array< OneD, const int > &	GetGlobalToUniversalMapUnique ()

NekDouble	GetLocalToGlobalSign (const int i) const

const Array< OneD, NekDouble > &	GetLocalToGlobalSign () const

void	LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

void	LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

void	GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const

void	GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const

void	Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

void	Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

void	UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const

void	UniversalAssemble (NekVector< NekDouble > &pGlobal) const

void	UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const

int	GetLocalToGlobalBndMap (const int i) const
	Retrieve the global index of a given local boundary mode. 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...

int	GetBndCondCoeffsToGlobalCoeffsMap (const int i)
	Retrieves the global index corresponding to a boundary expansion mode. More...

const Array< OneD, const int > &	GetBndCondCoeffsToGlobalCoeffsMap ()
	Retrieves the global indices corresponding to the boundary expansion modes. More...

NekDouble	GetBndCondCoeffsToGlobalCoeffsSign (const int i)
	Returns the modal sign associated with a given boundary expansion mode. More...

int	GetBndCondTraceToGlobalTraceMap (const int i)
	Returns the global index of the boundary trace giving the index on the boundary expansion. More...

const Array< OneD, const int > &	GetBndCondTraceToGlobalTraceMap ()

int	GetNumGlobalDirBndCoeffs () const
	Returns the number of global Dirichlet boundary coefficients. More...

int	GetNumLocalDirBndCoeffs () const
	Returns the number of local Dirichlet boundary coefficients. More...

int	GetNumGlobalBndCoeffs () const
	Returns the total number of global boundary coefficients. More...

int	GetNumLocalBndCoeffs () const
	Returns the total number of local boundary coefficients. More...

int	GetNumLocalCoeffs () const
	Returns the total number of local coefficients. More...

int	GetNumGlobalCoeffs () const
	Returns the total number of global coefficients. More...

bool	GetSingularSystem () const
	Retrieves if the system is singular (true) or not (false) More...

void	GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const

void	GlobalToLocalBnd (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const

void	GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc, int offset) const

void	GlobalToLocalBnd (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const

void	LocalBndToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const

void	LocalBndToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

void	LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const

void	LocalBndToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

void	AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const

void	AssembleBnd (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

void	AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset) const

void	AssembleBnd (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

void	UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal) const

void	UniversalAssembleBnd (NekVector< NekDouble > &pGlobal) const

void	UniversalAssembleBnd (Array< OneD, NekDouble > &pGlobal, int offset) const

int	GetFullSystemBandWidth () const

int	GetNumNonDirVertexModes () const

int	GetNumNonDirEdgeModes () const

int	GetNumNonDirFaceModes () const

int	GetNumDirEdges () const

int	GetNumDirFaces () const

int	GetNumNonDirEdges () const

int	GetNumNonDirFaces () const

void	PrintStats (std::ostream &out, std::string variable) const

const Array< OneD, const int > &	GetExtraDirEdges ()

boost::shared_ptr< AssemblyMap >	LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)

int	GetBndSystemBandWidth () const
	Returns the bandwidth of the boundary system. More...

int	GetStaticCondLevel () const
	Returns the level of static condensation for this map. More...

int	GetNumPatches () const
	Returns the number of patches in this static condensation level. More...

const Array< OneD, const unsigned int > &	GetNumLocalBndCoeffsPerPatch ()
	Returns the number of local boundary coefficients in each patch. More...

const Array< OneD, const unsigned int > &	GetNumLocalIntCoeffsPerPatch ()
	Returns the number of local interior coefficients in each patch. More...

const AssemblyMapSharedPtr	GetNextLevelLocalToGlobalMap () const
	Returns the local to global mapping for the next level in the multi-level static condensation. More...

void	SetNextLevelLocalToGlobalMap (AssemblyMapSharedPtr pNextLevelLocalToGlobalMap)

const PatchMapSharedPtr &	GetPatchMapFromPrevLevel (void) const
	Returns the patch map from the previous level of the multi-level static condensation. More...

bool	AtLastLevel () const
	Returns true if this is the last level in the multi-level static condensation. More...

GlobalSysSolnType	GetGlobalSysSolnType () const
	Returns the method of solving global systems. More...

PreconditionerType	GetPreconType () const

NekDouble	GetIterativeTolerance () const

int	GetSuccessiveRHS () const

int	GetLowestStaticCondLevel () const

Additional Inherited Members
Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMapCG
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, set< int > &extraDirVerts, set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)

void	SetUpUniversalC0ContMap (const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)

void	CalculateFullSystemBandWidth ()
	Calculate the bandwith of the full matrix system. More...

virtual int	v_GetLocalToGlobalMap (const int i) const

virtual int	v_GetGlobalToUniversalMap (const int i) const

virtual int	v_GetGlobalToUniversalMapUnique (const int i) const

virtual const Array< OneD, const int > &	v_GetLocalToGlobalMap ()

virtual const Array< OneD, const int > &	v_GetGlobalToUniversalMap ()

virtual const Array< OneD, const int > &	v_GetGlobalToUniversalMapUnique ()

virtual NekDouble	v_GetLocalToGlobalSign (const int i) const

virtual const Array< OneD, NekDouble > &	v_GetLocalToGlobalSign () const

virtual void	v_LocalToGlobal (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

virtual void	v_LocalToGlobal (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

virtual void	v_GlobalToLocal (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const

virtual void	v_GlobalToLocal (const NekVector< NekDouble > &global, NekVector< NekDouble > &loc) const

virtual void	v_Assemble (const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const

virtual void	v_Assemble (const NekVector< NekDouble > &loc, NekVector< NekDouble > &global) const

virtual void	v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal) const

virtual void	v_UniversalAssemble (NekVector< NekDouble > &pGlobal) const

virtual void	v_UniversalAssemble (Array< OneD, NekDouble > &pGlobal, int offset) const

virtual int	v_GetFullSystemBandWidth () const

virtual int	v_GetNumNonDirVertexModes () const

virtual int	v_GetNumNonDirEdgeModes () const

virtual int	v_GetNumNonDirFaceModes () const

virtual int	v_GetNumDirEdges () const

virtual int	v_GetNumDirFaces () const

virtual int	v_GetNumNonDirEdges () const

virtual int	v_GetNumNonDirFaces () const

virtual const Array< OneD, const int > &	v_GetExtraDirEdges ()

virtual AssemblyMapSharedPtr	v_LinearSpaceMap (const ExpList &locexp, GlobalSysSolnType solnType)
	Construct an AssemblyMapCG object which corresponds to the linear space of the current object. More...

Protected Member Functions inherited from Nektar::MultiRegions::AssemblyMap
void	CalculateBndSystemBandWidth ()
	Calculates the bandwidth of the boundary system. More...

void	GlobalToLocalBndWithoutSign (const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc)

Protected Attributes inherited from Nektar::MultiRegions::AssemblyMapCG
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...

map< int, vector< ExtraDirDof > >	m_extraDirDofs
	Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor. More...

Protected Attributes inherited from Nektar::MultiRegions::AssemblyMap
LibUtilities::SessionReaderSharedPtr	m_session
	Session object. More...

LibUtilities::CommSharedPtr	m_comm
	Communicator. More...

size_t	m_hash
	Hash for map. More...

int	m_numLocalBndCoeffs
	Number of local boundary coefficients. More...

int	m_numGlobalBndCoeffs
	Total number of global boundary coefficients. More...

int	m_numLocalDirBndCoeffs
	Number of Local Dirichlet Boundary Coefficients. More...

int	m_numGlobalDirBndCoeffs
	Number of Global Dirichlet Boundary Coefficients. More...

bool	m_systemSingular
	Flag indicating if the system is singular or not. More...

int	m_numLocalCoeffs
	Total number of local coefficients. More...

int	m_numGlobalCoeffs
	Total number of global coefficients. More...

bool	m_signChange
	Flag indicating if modes require sign reversal. More...

Array< OneD, int >	m_localToGlobalBndMap
	Integer map of local boundary coeffs to global space. More...

Array< OneD, NekDouble >	m_localToGlobalBndSign
	Integer sign of local boundary coeffs to global space. More...

Array< OneD, int >	m_bndCondCoeffsToGlobalCoeffsMap
	Integer map of bnd cond coeffs to global coefficients. More...

Array< OneD, NekDouble >	m_bndCondCoeffsToGlobalCoeffsSign
	Integer map of bnd cond coeffs to global coefficients. More...

Array< OneD, int >	m_bndCondTraceToGlobalTraceMap
	Integer map of bnd cond trace number to global trace number. More...

Array< OneD, int >	m_globalToUniversalBndMap
	Integer map of process coeffs to universal space. More...

Array< OneD, int >	m_globalToUniversalBndMapUnique
	Integer map of unique process coeffs to universal space (signed) More...

GlobalSysSolnType	m_solnType
	The solution type of the global system. More...

int	m_bndSystemBandWidth
	The bandwith of the global bnd system. More...

PreconditionerType	m_preconType
	Type type of preconditioner to use in iterative solver. More...

NekDouble	m_iterativeTolerance
	Tolerance for iterative solver. More...

int	m_successiveRHS
	sucessive RHS for iterative solver More...

Gs::gs_data *	m_gsh

Gs::gs_data *	m_bndGsh

int	m_staticCondLevel
	The level of recursion in the case of multi-level static condensation. More...

int	m_numPatches
	The number of patches (~elements) in the current level. More...

Array< OneD, unsigned int >	m_numLocalBndCoeffsPerPatch
	The number of bnd dofs per patch. More...

Array< OneD, unsigned int >	m_numLocalIntCoeffsPerPatch
	The number of int dofs per patch. More...

AssemblyMapSharedPtr	m_nextLevelLocalToGlobalMap
	Map from the patches of the previous level to the patches of the current level. More...

int	m_lowestStaticCondLevel
	Lowest static condensation level. More...

Detailed Description

Definition at line 45 of file CoupledLocalToGlobalC0ContMap.h.

Constructor & Destructor Documentation

Nektar::CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const SpatialDomains::MeshGraphSharedPtr &	graph,
		const SpatialDomains::BoundaryConditionsSharedPtr &	boundaryConditions,
		const Array< OneD, MultiRegions::ExpListSharedPtr > &	fields,
		const MultiRegions::ExpListSharedPtr &	pressure,
		const int	nz_loc,
		const bool	CheckforSingularSys = `true`
	)

This is an vector extension of MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the Linearised Navier Stokes problem

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

: Fix this so that we can extract normals from edges

STEP 2a: Set the mean pressure modes to edges depending on type of direct solver technique;

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

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

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

STEP 4: Now, all ingredients are ready to set up the actual local to global mapping.

The remainder of the map consists of the element-interior degrees of freedom. This leads to the block-diagonal submatrix as each element-interior mode is globally orthogonal to modes in all other elements.

STEP 5: The boundary condition mapping is generated from the same vertex renumbering and fill in a unique interior map.

Definition at line 51 of file CoupledLocalToGlobalC0ContMap.cpp.

                                        :
         AssemblyMapCG(pSession)
     {
         int i,j,k,n;
         int cnt = 0,offset=0;
         int meshVertId;
         int meshEdgeId;
         int bndEdgeCnt;
         int globalId;
         int nEdgeCoeffs;
         int nEdgeInteriorCoeffs;
         int firstNonDirGraphVertId;
         int nLocBndCondDofs = 0;
         int nLocDirBndCondDofs = 0;
         int nExtraDirichlet = 0;
         StdRegions::StdExpansion2DSharedPtr locExpansion;
         LocalRegions::SegExpSharedPtr       bndSegExp;
         LibUtilities::BasisType             bType;
         StdRegions::Orientation         edgeOrient;
         Array<OneD, unsigned int>           edgeInteriorMap;
         Array<OneD, int>                    edgeInteriorSign;
         int nvel = fields.num_elements();
         
         const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
         int eid, id, diff;
         int nel = fields[0]->GetNumElmts();
 
         MultiRegions::PeriodicMap periodicVerts;
         MultiRegions::PeriodicMap periodicEdges;
         MultiRegions::PeriodicMap periodicFaces;
         vector<map<int,int> > ReorderedGraphVertId(3);
         MultiRegions::BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
         int staticCondLevel = 0;
 
         if(CheckforSingularSys) //all singularity checking by setting flag to true
         {
             m_systemSingular = true;
         }
         else  // Turn off singular checking by setting flag to false
         {
             m_systemSingular = false;
         }
         
         /**
          * STEP 1: Wrap boundary conditions vector in an array
          * (since routine is set up for multiple fields) and call
          * the graph re-odering subroutine to obtain the reordered
          * values
          */
 
         // Obtain any periodic information and allocate default mapping array
         fields[0]->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
 
 
         const Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp = fields[0]->GetBndCondExpansions();
 
         Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> > bndConditionsVec(nvel);
         for(i = 0; i < nvel; ++i)
         {
             bndConditionsVec[i] = fields[i]->GetBndConditions(); 
         }
         
         map<int,int> IsDirVertDof;
         map<int,int> IsDirEdgeDof;
         map<int,int>::iterator mapIt;
 
         SpatialDomains::Geometry1DSharedPtr g;
         for(j = 0; j < bndCondExp.num_elements(); ++j)
         {
             map<int,int> BndExpVids;
             // collect unique list of vertex ids for this expansion
             for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
             {
                 g = bndCondExp[j]->GetExp(k)->as<LocalRegions::Expansion1D>()
                                             ->GetGeom1D();
                 BndExpVids[g->GetVid(0)] = g->GetVid(0);
                 BndExpVids[g->GetVid(1)] = g->GetVid(1);
             }
             
             for(i = 0; i < nvel; ++i)
             {
                 if(bndConditionsVec[i][j]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
                 {
                     // set number of Dirichlet conditions along edge
                     for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
                     {
                         IsDirEdgeDof[bndCondExp[j]->GetExp(k)
                                      ->as<LocalRegions::Expansion1D>()
                                      ->GetGeom1D()->GetEid()] += 1;
                     }
                     
                             
                     // Set number of Dirichlet conditions at vertices
                     // with a clamp on its maximum value being nvel to
                     // handle corners between expansions
                     for(mapIt = BndExpVids.begin(); mapIt !=  BndExpVids.end(); mapIt++)
                     {
                         id = IsDirVertDof[mapIt->second]+1; 
                         IsDirVertDof[mapIt->second] = (id > nvel)?nvel:id;
                     }
                 }
                 else
                 { 
                     // Check to see that edge normals have non-zero
                     // component in this direction since otherwise
                     // also can be singular.
                     /// @TODO: Fix this so that we can extract normals from edges
                     for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
                     {
                         Array<OneD,Array<OneD,NekDouble> > locnorm;
                         LocalRegions::Expansion1DSharedPtr loc_exp
                             = bndCondExp[j]->GetExp(k)
                                 ->as<LocalRegions::Expansion1D>();
                         locnorm = loc_exp->GetLeftAdjacentElementExp()->GetEdgeNormal(loc_exp->GetLeftAdjacentElementEdge());
                         //locnorm = bndCondExp[j]->GetExp(k)->Get GetMetricInfo()->GetNormal();
                         
                         int ndir = locnorm.num_elements();
                         if(i < ndir) // account for Fourier version where n can be larger then ndir
                         {
                             for(int l = 0; l < locnorm[0].num_elements(); ++l)
                             {
                                 if(fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
                                 {
                                     m_systemSingular = false;
                                     break;
                                 }
                             }
                         }
                         if(m_systemSingular == false)
                         {
                             break;
                         }
                     }
                 }
             }
         }
 
         Array<OneD, map<int,int> >Dofs(2);
 
         Array<OneD, int> AddMeanPressureToEdgeId(nel,-1);
         int edgeId,vertId;
 
 
         // special case of singular problem - need to fix one pressure
         // dof to a dirichlet edge. Since we attached pressure dof to
         // last velocity component of edge need to make sure this
         // component is Dirichlet
         if(m_systemSingular)
         {
             id  = -1;
             for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
             {
                 if(bndConditionsVec[nvel-1][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
                 {
                     id = bndCondExp[i]->GetExp(0)
                                 ->as<LocalRegions::Expansion1D>()->GetGeom1D()
                                 ->GetEid();
                     break;
                 }
             }
             
             ASSERTL0(id != -1," Did not find an edge to attach singular pressure degree of freedom");
 
             // determine element with this edge id. There may be a
             // more direct way of getting element from spatialDomains
             for(i = 0; i < nel; ++i)
             {
                 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
                 {
                     edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
                                              ->GetGeom2D())->GetEid(j);
                     
                     if(edgeId == id)
                     {
                         AddMeanPressureToEdgeId[i] = id;
                         break;
                     }
                 }
                 
                 if(AddMeanPressureToEdgeId[i] != -1)
                 {
                     break;
                 }
             }
         }
 
 
         for(i = 0; i < nel; ++i)
         {
             eid = fields[0]->GetOffset_Elmt_Id(i);
             for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
             {
                 vertId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
                                            ->GetGeom2D())->GetVid(j);
                 if(Dofs[0].count(vertId) == 0)
                 {
                     Dofs[0][vertId] = nvel*nz_loc;
                     
                     // Adjust for a Dirichlet boundary condition to give number to be solved
                     if(IsDirVertDof.count(vertId) != 0)
                     {
                         Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc; 
                     }
                 }
 
                 edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
                                            ->GetGeom2D())->GetEid(j);
                 if(Dofs[1].count(edgeId) == 0)
                 {
                     Dofs[1][edgeId] = nvel*(locExpVector[eid]->GetEdgeNcoeffs(j)-2)*nz_loc;
                 }
                 
                 // Adjust for Dirichlet boundary conditions to give number to be solved
                 if(IsDirEdgeDof.count(edgeId) != 0)
                 {
                     Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[eid]->GetEdgeNcoeffs(j)-2);
                 }
             }
         }
 
         set<int> extraDirVerts, extraDirEdges;
 
         CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false,
                     periodicVerts, periodicEdges, periodicFaces,
                     ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
                     extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
         /*
         SetUp2DGraphC0ContMap(*fields[0],
                               bndCondExp,
                               bndConditionsVec,
                               periodicVerts,          periodicEdges,
                               Dofs,                   ReorderedGraphVertId,
                               firstNonDirGraphVertId, nExtraDirichlet,
                               bottomUpGraph, extraDir,  false,  4);
         */
         
         /**
          * STEP 2a: Set the mean pressure modes to edges depending on
          * type of direct solver technique;
          */
         
         // determine which edge to add mean pressure dof based on
         // ensuring that at least one pressure dof from an internal
         // patch is associated with its boundary system
         if(m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
         {
 
             
             FindEdgeIdToAddMeanPressure(ReorderedGraphVertId,
                  nel, locExpVector,
                  edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
                  bottomUpGraph,
                  AddMeanPressureToEdgeId);
         }   
         
         // Set unset elmts to non-Dirichlet edges. 
         // special case of singular problem - need to fix one
         // pressure dof to a dirichlet edge 
         for(i = 0; i < nel; ++i)
         {
             eid = fields[0]->GetOffset_Elmt_Id(i);
             for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
             {
                 edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
                                            ->GetGeom2D())->GetEid(j);
 
                 if(IsDirEdgeDof.count(edgeId) == 0) // interior edge
                 {
                     // setup AddMeanPressureToEdgeId to decide where to
                     // put pressure
                     if(AddMeanPressureToEdgeId[eid] == -1)
                     {
                         AddMeanPressureToEdgeId[eid] = edgeId;
                     }
                 }
             }
             ASSERTL0((AddMeanPressureToEdgeId[eid] != -1),"Did not determine "
                      "an edge to attach mean pressure dof");
             // Add the mean pressure degree of freedom to this edge
             Dofs[1][AddMeanPressureToEdgeId[eid]] += nz_loc;
         }
         
         map<int,int> pressureEdgeOffset;
                 
         /**
          * STEP 2: Count out the number of Dirichlet vertices and edges first
          */
         for(i = 0; i < bndCondExp.num_elements(); i++)
         {
             for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
             {
                 bndSegExp = bndCondExp[i]->GetExp(j)
                                          ->as<LocalRegions::SegExp>();
                 for(k = 0; k < nvel; ++k)
                 {
                     if(bndConditionsVec[k][i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
                     {
                         nLocDirBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
                     }
                     nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
                 }
             }
         }
 
         if(m_systemSingular)
         {
             m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet+nz_loc;
         }
         else
         {
             m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet;
         }
 
         /**
          * STEP 3: Set up an array which contains the offset information of
          * the different graph vertices.
          *
          * This basically means to identify how many global degrees of
          * freedom the individual graph vertices correspond. Obviously,
          * the graph vertices corresponding to the mesh-vertices account
          * for a single global DOF. However, the graph vertices
          * corresponding to the element edges correspond to 2*(N-2) global DOF
          * where N is equal to the number of boundary modes on this edge.
          */
         Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
         graphVertOffset[0] = 0;
         
         m_signChange = false;
 
         for(i = 0; i < nel; ++i)
         {
             eid = fields[0]->GetOffset_Elmt_Id(i);
             locExpansion = locExpVector[eid]->as<StdRegions::StdExpansion2D>();
 
             for(j = 0; j < locExpansion->GetNedges(); ++j)
             {
                 nEdgeCoeffs = locExpansion->GetEdgeNcoeffs(j);
                 meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
                                           ->GetGeom2D())->GetEid(j);
                 meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
                                           ->GetGeom2D())->GetVid(j);
 
                 for(k = 0; k < nvel*nz_loc; ++k)
                 {
                     graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+k] = 1;
                     graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+k] = (nEdgeCoeffs-2);
                 }
                 
                 bType = locExpansion->GetEdgeBasisType(j);
                 // need a sign vector for modal expansions if nEdgeCoeffs >=4
                 if( (nEdgeCoeffs >= 4)&&
                     ( (bType == LibUtilities::eModified_A)||
                       (bType == LibUtilities::eModified_B) ) )
                 {
                     m_signChange = true;
                 }
             }
         }
 
         // Add mean pressure modes; 
         for(i = 0; i < nel; ++i)
         {
             graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1] += nz_loc;
             //graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc] += nz_loc;
         }
         
         // Negate the vertices and edges with only a partial
         // Dirichlet conditon. Essentially we check to see if an edge
         // has a mixed Dirichlet with Neumann/Robin Condition and if
         // so negate the offset associated with this vertex. 
         
         map<int,int> DirVertChk;
 
         for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
         {
             cnt = 0;
             for(j = 0; j < nvel; ++j)
             {
                 if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
                 {
                     cnt ++;
                 }
             }
             
             // Case where partial Dirichlet boundary condition
             if((cnt > 0)&&(cnt < nvel)) 
             {
                 for(j  = 0; j < nvel; ++j)
                 {
                     if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
                     {
                         //negate graph offsets which should be
                         //Dirichlet conditions
                         for(k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
                         {
                             // vertices with mix condition;
                             id = bndCondExp[i]->GetExp(k)
                                               ->as<LocalRegions::Expansion1D>()
                                               ->GetGeom1D()->GetVid(0);
                             if(DirVertChk.count(id*nvel+j) == 0)
                             {
                                 DirVertChk[id*nvel+j] = 1;
                                 for(n = 0; n < nz_loc; ++n)
                                 {
                                     graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc + n] *= -1; 
                                 }
                             }
 
                             id = bndCondExp[i]->GetExp(k)
                                               ->as<LocalRegions::Expansion1D>()
                                               ->GetGeom1D()->GetVid(1);
                             if(DirVertChk.count(id*nvel+j) == 0)
                             {
                                 DirVertChk[id*nvel+j] = 1;
                                 for(n = 0; n < nz_loc; ++n)
                                 {
                                     graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc+n] *= -1;
                                 }
                             }
                             
                             // edges with mixed id; 
                             id = bndCondExp[i]->GetExp(k)
                                               ->as<LocalRegions::Expansion1D>()
                                               ->GetGeom1D()->GetEid();
                             for(n = 0; n < nz_loc; ++n)
                             {
                                 graphVertOffset[ReorderedGraphVertId[1][id]*nvel*nz_loc+j*nz_loc +n] *= -1; 
                             }
                         }
                     }
                 }
             }
         }
         
 
         cnt = 0;
         // assemble accumulative list of full Dirichlet values. 
         for(i = 0; i < firstNonDirGraphVertId*nvel*nz_loc; ++i)
         {
             diff = abs(graphVertOffset[i]); 
             graphVertOffset[i] = cnt; 
             cnt += diff;
         }
         
         // set Dirichlet values with negative values to Dirichlet value
         for(i = firstNonDirGraphVertId*nvel*nz_loc; i <  graphVertOffset.num_elements(); ++i)
         {
             if(graphVertOffset[i] < 0)
             {
                 diff = -graphVertOffset[i]; 
                 graphVertOffset[i] = -cnt; 
                 cnt += diff;
             }
         }
     
         // Accumulate all interior degrees of freedom with positive values
         m_numGlobalDirBndCoeffs = cnt;
 
         // offset values
         for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
         {
             if(graphVertOffset[i] >= 0)
             {
                 diff = graphVertOffset[i]; 
                 graphVertOffset[i] = cnt; 
                 cnt += diff; 
             }
         }
 
         // Finally set negative entries (corresponding to Dirichlet
         // values ) to be positive
         for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
         {
             if(graphVertOffset[i] < 0)
             {
                 graphVertOffset[i] = -graphVertOffset[i]; 
             }
         }
 
 
         // Allocate the proper amount of space for the class-data and fill
         // information that is already known
         cnt = 0; 
         m_numLocalBndCoeffs = 0;
         m_numLocalCoeffs = 0;
 
         for(i = 0; i < nel; ++i)
         {
             m_numLocalBndCoeffs += nz_loc*(nvel*locExpVector[i]->NumBndryCoeffs() + 1);
             // add these coeffs up separately since
             // pressure->GetNcoeffs can include the coefficient in
             // multiple planes.
             m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs()-1)*nz_loc;
         }
 
         m_numLocalCoeffs += m_numLocalBndCoeffs; 
 
 
         m_localToGlobalMap    = Array<OneD, int>(m_numLocalCoeffs,-1);
         m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
         m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(nLocBndCondDofs,-1);
         
 
         // Set default sign array. 
         m_localToGlobalSign    = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
         m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
 
         m_staticCondLevel = staticCondLevel;
         m_numPatches = nel;
         
         m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
         m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
 
         for(i = 0; i < nel; ++i)
         {
             m_numLocalBndCoeffsPerPatch[i] = (unsigned int) nz_loc*(nvel*locExpVector[fields[0]->GetOffset_Elmt_Id(i)]->NumBndryCoeffs() + 1);
             m_numLocalIntCoeffsPerPatch[i] = (unsigned int) nz_loc*(pressure->GetExp(i)->GetNcoeffs()-1);
         }
         
         /**
          * STEP 4: Now, all ingredients are ready to set up the actual
          * local to global mapping.
          *
          * The remainder of the map consists of the element-interior
          * degrees of freedom. This leads to the block-diagonal submatrix
          * as each element-interior mode is globally orthogonal to modes
          * in all other elements.
          */
         cnt = 0;
         int nv,velnbndry;
         Array<OneD, unsigned int> bmap;
         
        
         // Loop over all the elements in the domain in shuffled
         // ordering (element type consistency)
         for(i = 0; i < nel; ++i)
         {
             eid = fields[0]->GetOffset_Elmt_Id(i);
             locExpansion = locExpVector[eid]->as<StdRegions::StdExpansion2D>();
 
             velnbndry = locExpansion->NumBndryCoeffs();
 
             // require an inverse ordering of the bmap system to store
             // local numbering system which takes matrix these
             // matrices. Therefore get hold of elemental bmap and set
             // up an inverse map
             map<int,int> inv_bmap;
             locExpansion->GetBoundaryMap(bmap);
             for(j = 0; j < bmap.num_elements(); ++j)
             {
                 inv_bmap[bmap[j]] = j;
             }
             
             // Loop over all edges (and vertices) of element i
             for(j = 0; j < locExpansion->GetNedges(); ++j)
             {
                 nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
                 edgeOrient   = (locExpansion->as<LocalRegions::Expansion2D>()
                                             ->GetGeom2D())->GetEorient(j);
                 meshEdgeId   = (locExpansion->as<LocalRegions::Expansion2D>()
                                             ->GetGeom2D())->GetEid(j);
                 meshVertId   = (locExpansion->as<LocalRegions::Expansion2D>()
                                             ->GetGeom2D())->GetVid(j);
                 
                 locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
                 // Set the global DOF for vertex j of element i
                 
                 for(nv = 0; nv < nvel*nz_loc; ++nv)
                 {
                     m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
                     // Set the global DOF's for the interior modes of edge j
                     for(k = 0; k < nEdgeInteriorCoeffs; ++k)
                     {
                         m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] =  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
                     }
                 }
                 
                 // Fill the sign vector if required
                 if(m_signChange)
                 {
                     for(nv = 0; nv < nvel*nz_loc; ++nv)
                     {
                         for(k = 0; k < nEdgeInteriorCoeffs; ++k)
                         {
                             m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (NekDouble) edgeInteriorSign[k];
                         }
                     }
                 }
             }
 
             // use difference between two edges of the AddMeanPressureEdgeId to det nEdgeInteriorCoeffs. 
             nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc];
 
             int psize = pressure->GetExp(eid)->GetNcoeffs();
             for(n = 0; n < nz_loc; ++n)
             {
                 m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[eid]];
                 
                 pressureEdgeOffset[AddMeanPressureToEdgeId[eid]] += 1;
             }
             
             cnt += (velnbndry*nvel+ psize)*nz_loc;
         }
        
         // Set up the mapping for the boundary conditions
         offset = cnt = 0;
         for(nv = 0; nv < nvel; ++nv)
         {
             for(i = 0; i < bndCondExp.num_elements(); i++)
             {
                 for(n = 0; n < nz_loc; ++n)
                 {
                     int ncoeffcnt = 0;
                     for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
                     {
                         bndSegExp  = bndCondExp[i]->GetExp(j)
                                                   ->as<LocalRegions::SegExp>();
 
                         cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
                         for(k = 0; k < 2; k++)
                         {
                             meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
                             m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
                         }
                     
                         meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
                         bndEdgeCnt = 0;
                         nEdgeCoeffs = bndSegExp->GetNcoeffs();
                         for(k = 0; k < nEdgeCoeffs; k++)
                         {
                             if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
                             {
                                 m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
                                     graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
                                 bndEdgeCnt++;
                             }
                         }
                         ncoeffcnt += nEdgeCoeffs; 
                     }
                     // Note: Can not use bndCondExp[i]->GetNcoeffs()
                     // due to homogeneous extension not returning just
                     // the value per plane
                     offset += ncoeffcnt; 
                 }
             }
         }
         
         globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
         m_numGlobalBndCoeffs = globalId; 
         
         /**
          * STEP 5: The boundary condition mapping is generated from the
          * same vertex renumbering and fill in a unique interior map.
          */
         cnt=0;
         for(i = 0; i < m_numLocalCoeffs; ++i)
         {
             if(m_localToGlobalMap[i] == -1)
             {
                 m_localToGlobalMap[i] = globalId++;
             }
             else
             {
                 if(m_signChange)
                 {
                     m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
                 }
                 m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
             }
         }
         m_numGlobalCoeffs = globalId;
 
         // Set up the local to global map for the next level when using
         // multi-level static condensation
         if( m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond) )
         {
             if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
             {
                 Array<OneD, int> vwgts_perm(
                     Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
                 for(i = 0; i < locExpVector.size(); ++i)
                 {
                     eid = fields[0]->GetOffset_Elmt_Id(i);
                     locExpansion = locExpVector[eid]
                                             ->as<StdRegions::StdExpansion2D>();
                     for(j = 0; j < locExpansion->GetNverts(); ++j)
                     {
                         meshEdgeId = (locExpansion
                                         ->as<LocalRegions::Expansion2D>()
                                         ->GetGeom2D())->GetEid(j);
                         meshVertId = (locExpansion
                                         ->as<LocalRegions::Expansion2D>()
                                         ->GetGeom2D())->GetVid(j);
                             
                         if(ReorderedGraphVertId[0][meshVertId] >= 
                            firstNonDirGraphVertId)
                         {
                             vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
                                        firstNonDirGraphVertId] = 
                                 Dofs[0][meshVertId];
                         }
                             
                         if(ReorderedGraphVertId[1][meshEdgeId] >= 
                            firstNonDirGraphVertId)
                         {
                             vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
                                        firstNonDirGraphVertId] = 
                                 Dofs[1][meshEdgeId];
                         }
                     }
                 }
 
                 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
                     
                 m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
                     AllocateSharedPtr(this,bottomUpGraph);
             }
         }
     }

Member Function Documentation

void CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure	(	vector< map< int, int > > &	ReorderedGraphVertId,
		int &	nel,
		const LocalRegions::ExpansionVector &	locExpVector,
		int &	edgeId,
		int &	vertId,
		int &	firstNonDirGraphVertId,
		map< int, int > &	IsDirEdgeDof,
		MultiRegions::BottomUpSubStructuredGraphSharedPtr &	bottomUpGraph,
		Array< OneD, int > &	AddMeanPressureToEdgeId
	)

Definition at line 780 of file CoupledLocalToGlobalC0ContMap.cpp.

References Nektar::iterator.

Referenced by CoupledLocalToGlobalC0ContMap().

 {
     
     int i,j,k;
         
     // Make list of homogeneous graph edges to elmt mappings
     Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(),2,-1);
     map<int,int> HomGraphEdgeIdToEdgeId;
     
     for(i = 0; i < nel; ++i)
     {
         for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
         {
             edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
                                      ->GetGeom2D())->GetEid(j);
             
             // note second condition stops us using mixed boundary condition 
             if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
                && (IsDirEdgeDof.count(edgeId) == 0))
             {
                 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
                 
                 if(EdgeIdToElmts[edgeId][0] == -1)
                 {
                     EdgeIdToElmts[edgeId][0] = i;
                 }
                 else
                 {
                     EdgeIdToElmts[edgeId][1] = i;
                 }
             }
         }
     }
     
     
     map<int,int>::iterator mapIt;
     
     // Start at second to last level and find edge on boundary
     // to attach element
     int nlevels = bottomUpGraph->GetNlevels();
     
     // determine a default edge to attach pressure modes to
     // which is part of the inner solve;
     int defedge = -1;
     
     vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
     for(i = 0; i < bndgraphs.size(); ++i)
     {
         int GlobIdOffset = bndgraphs[i]->GetIdOffset();
         
         for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
         {
             // find edge in graph vert list
             if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
             {
                 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
 
                 if(defedge == -1)
                 {
                     defedge = edgeId;
                     break;
                 }
             }
         }
         if(defedge != -1)
         {
             break;
         }
     }
     
     for(int n = 1; n < nlevels; ++n)
     {
         // produce a map with a key that is the element id
         // that contains which next level patch it belongs to
         vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
         
         // Fill next level graph  of adjacent elements and their level
         map<int,int> ElmtInBndry;
         
         for(i = 0; i < bndgraphs.size(); ++i)
         {
             int GlobIdOffset = bndgraphs[i]->GetIdOffset();
             
             for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
             {
                 // find edge in graph vert list
                 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
                 {
                     edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];  
                     
                     if(EdgeIdToElmts[edgeId][0] != -1)
                     {
                         ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
                     }
                     if(EdgeIdToElmts[edgeId][1] != -1)
                     {
                         ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
                     }
                 }
             }
         }
         
         // Now search interior patches in this level for edges
         // that share the same element as a boundary edge and
         // assign this elmt that boundary edge
         vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
         for(i = 0; i < intgraphs.size(); ++i)
         {
             int GlobIdOffset = intgraphs[i]->GetIdOffset();
             bool SetEdge = false; 
             int elmtid = 0;
             for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
             {
                 // Check to see if graph vert is an edge 
                 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
                 {
                     edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
                     
                     for(k = 0; k < 2; ++k)
                     {
                         // relevant edge id
                         elmtid = EdgeIdToElmts[edgeId][k];
                         
                         if(elmtid != -1)
                         {
                             mapIt = ElmtInBndry.find(elmtid);
                             
                             if(mapIt != ElmtInBndry.end())
                             {
                                 // now find a edge in the next level boundary graph 
                                 int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
                                 for(int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
                                 {
                                     // find edge in graph vert list
                                     if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
                                     {
                                         //June 2012: commenting this condition apparently 
                                         //solved the bug caused by the edge reordering procedure
                                                                             
                                         //if(AddMeanPressureToEdgeId[elmtid] == -1)
                                         //{
                                         
                                         //AddMeanPressureToEdgeId[elmtid] = HomGraphEdgeIdToEdgeId[GlobIdOffset1+l];
                                         AddMeanPressureToEdgeId[elmtid] = defedge;
                                         
                                         //}
                                         SetEdge = true;
                                         break;
                                     }
                                 }
                             }
                         }
                     }
                 }
             }
             
             
             // if we have failed to find matching edge in next
             // level patch boundary then set last found elmt
             // associated to this interior patch to the
             // default edget value
             if(SetEdge == false)
             {
                 if(elmtid == -1) // find an elmtid in patch 
                 {
                     for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
                     {
                         if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
                         {
                             edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
                             for(k = 0; k < 2; ++k)
                             {
                                 // relevant edge id
                                 elmtid = EdgeIdToElmts[edgeId][k];
                                 if(elmtid != -1)
                                 {
                                     break;
                                 }
                             }
                         }
                         if(elmtid != -1)
                         {
                             break;
                         }
                     }
                 }
                 if(AddMeanPressureToEdgeId[elmtid] == -1)
                 {
                     AddMeanPressureToEdgeId[elmtid] = defedge;
                 }
             }
         }
     }
     
 }

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation