49 #include <boost/config.hpp>
50 #include <boost/graph/adjacency_list.hpp>
51 #include <boost/graph/cuthill_mckee_ordering.hpp>
52 #include <boost/graph/properties.hpp>
53 #include <boost/graph/bandwidth.hpp>
59 namespace MultiRegions
61 AssemblyMapDG::AssemblyMapDG():
62 m_numDirichletBndPhys(0)
78 const std::string variable):
81 int i, j, k, cnt, eid, id, id1, gid;
83 int nTraceExp = trace->GetExpSize();
84 int nbnd = bndCondExp.num_elements();
91 int nel = expList.size();
93 map<int, int> meshTraceId;
98 for(i = 0; i < nTraceExp; ++i)
100 meshTraceId[trace->GetExp(i)->GetGeom()->GetGlobalID()] = i;
105 for(i = 0; i < nel; ++i)
107 cnt += expList[i]->GetNtrace();
115 for(cnt = i = 0; i < nel; ++i)
119 for(j = 0; j < expList[i]->GetNtrace(); ++j)
121 id = expList[i]->GetGeom()->GetTid(j);
123 if(meshTraceId.count(
id) > 0)
126 trace->GetExp(meshTraceId.find(
id)->second);
130 ASSERTL0(
false,
"Failed to find trace map");
134 cnt += expList[i]->GetNtrace();
139 for(i = 0; i < nbnd; ++i)
141 cnt += bndCondExp[i]->GetExpSize();
150 for(i = 0; i < bndCondExp.num_elements(); ++i)
152 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
154 bndExp = bndCondExp[i]->GetExp(j);
155 traceGeom = bndExp->GetGeom();
156 id = traceGeom->GetGlobalID();
158 if(bndCond[i]->GetBoundaryConditionType() ==
179 for(i = 0; i < nel; ++i)
182 nbndry += expList[eid]->NumDGBndryCoeffs();
183 m_numLocalIntCoeffsPerPatch[i] = 0;
185 (
unsigned int) expList[eid]->NumDGBndryCoeffs();
201 typedef boost::adjacency_list<
202 boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
204 BoostGraph boostGraphObj;
205 int trace_id, trace_id1;
210 for(i = 0; i < nTraceExp; ++i)
212 id = trace->GetExp(i)->GetGeom()->GetGlobalID();
214 if (dirTrace.count(
id) == 0)
218 boost::add_vertex(boostGraphObj);
219 traceElmtGid[i] = cnt++;
224 traceElmtGid[i] = dirOffset;
225 dirOffset += trace->GetExp(i)->GetNcoeffs();
231 for(i = 0; i < nel; ++i)
235 for(j = 0; j < expList[eid]->GetNtrace(); ++j)
239 id = traceGeom->GetGlobalID();
240 trace_id = meshTraceId.find(
id)->second;
242 if(dirTrace.count(
id) == 0)
244 for(k = j+1; k < expList[eid]->GetNtrace(); ++k)
247 id1 = traceGeom->GetGlobalID();
248 trace_id1 = meshTraceId.find(id1)->second;
250 if(dirTrace.count(id1) == 0)
252 boost::add_edge((
size_t)traceElmtGid[trace_id],
253 (
size_t)traceElmtGid[trace_id1],
261 int nGraphVerts = nTraceExp - nDir;
267 for(i = 0; i < nGraphVerts; ++i)
269 vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
303 ASSERTL0(
false,
"Unrecognised solution type");
311 for(i = 0; i < nTraceExp - nDir; ++i)
313 traceElmtGid[perm[i]+nDir] = cnt;
314 cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
320 for(i = 0; i < nel; ++i)
326 for(j = 0; j < exp->GetNtrace(); ++j)
329 id = traceGeom->GetGlobalID();
330 gid = traceElmtGid[meshTraceId.find(
id)->second];
332 const int nDim = expList[eid]->GetNumBases();
341 order_e = expList[eid]->GetEdgeNcoeffs(j);
345 for(k = 0; k < order_e; ++k)
359 for (k = 2; k < order_e; ++k)
365 for(k = 3; k < order_e; k+=2)
367 m_localToGlobalBndSign[cnt+k] = -1.0;
374 for(k = 0; k < order_e; ++k)
383 for(k = 0; k < order_e; ++k)
391 ASSERTL0(
false,
"Boundary type not permitted");
398 order_e = expList[eid]->GetFaceNcoeffs(j);
400 std::map<int, int> orientMap;
411 expList[eid]->GetFaceToElementMap(j,fo,elmMap1,elmSign1);
412 expList[eid]->GetFaceToElementMap(
415 for (k = 0; k < elmMap1.num_elements(); ++k)
420 for (
int l = 0; l < elmMap2.num_elements(); ++l)
422 if (elmMap1[k] == elmMap2[l])
429 ASSERTL2(idx != -1,
"Problem with face to element map!");
433 for(k = 0; k < order_e; ++k)
436 m_localToGlobalBndSign[k+cnt] = elmSign2[orientMap[k]];
446 for(i = 0; i < nbnd; ++i)
448 cnt += bndCondExp[i]->GetNcoeffs();
454 int nbndexp = 0, bndOffset, bndTotal = 0;
455 for(cnt = i = 0; i < nbnd; ++i)
457 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
459 bndExp = bndCondExp[i]->GetExp(j);
460 id = bndExp->GetGeom()->GetGlobalID();
461 gid = traceElmtGid[meshTraceId.find(
id)->second];
462 bndOffset = bndCondExp[i]->GetCoeff_Offset(j) + bndTotal;
467 for(k = 0; k < bndExp->GetNcoeffs(); ++k)
473 nbndexp += bndCondExp[i]->GetExpSize();
474 bndTotal += bndCondExp[i]->GetNcoeffs();
492 for(
int i = 0; i < nGraphVerts; i++)
494 vwgts_perm[i] = vwgts[perm[i]];
497 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
505 for(i = 0; i < bndCondExp.num_elements(); ++i)
507 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
509 bndExp = bndCondExp[i]->GetExp(j);
510 traceGeom = bndExp->GetGeom();
511 id = traceGeom->GetGlobalID();
513 meshTraceId.find(
id)->second;
555 for(i = 0; i < locExpVector.size(); ++i)
557 locExpansion = locExpVector[i];
558 nDim = locExpansion->GetShapeDimension();
563 maxDof = (1 > maxDof ? 1 : maxDof);
567 for (j = 0; j < locExpansion->GetNedges(); ++j)
569 dof = locExpansion->GetEdgeNcoeffs(j);
570 maxDof = (dof > maxDof ? dof : maxDof);
575 for (j = 0; j < locExpansion->GetNfaces(); ++j)
577 dof = locExpansion->GetFaceNcoeffs(j);
578 maxDof = (dof > maxDof ? dof : maxDof);
586 for(i = 0; i < locExpVector.size(); ++i)
588 locExpansion = locExpVector[i];
589 nDim = locExpansion->GetShapeDimension();
598 int nverts = locExpansion->GetNverts();
599 for(j = 0; j < nverts; ++j)
603 id = locPointExp->
GetGeom()->GetGlobalID();
606 =
id * maxDof + j + 1;
612 for(j = 0; j < locExpansion->GetNedges(); ++j)
618 order_e = locExpVector[eid]->GetEdgeNcoeffs(j);
620 map<int,int> orientMap;
625 locExpVector[eid]->GetEdgeToElementMap(j, locExpVector[eid]->GetEorient(j), map2, sign2);
627 for (k = 0; k < map1.num_elements(); ++k)
632 for (
int l = 0; l < map2.num_elements(); ++l)
634 if (map1[k] == map2[l])
641 ASSERTL2(idx != -1,
"Problem with face to element map!");
645 for(k = 0; k < order_e; ++k)
649 =
id * maxDof + orientMap[k] + 1;
656 for(j = 0; j < locExpansion->GetNfaces(); ++j)
663 order_e = locExpVector[eid]->GetFaceNcoeffs(j);
665 map<int,int> orientMap;
670 locExpVector[eid]->GetFaceToElementMap(j, locExpVector[eid]->GetForient(j), map2, sign2);
672 for (k = 0; k < map1.num_elements(); ++k)
677 for (
int l = 0; l < map2.num_elements(); ++l)
679 if (map1[k] == map2[l])
686 ASSERTL2(idx != -1,
"Problem with face to element map!");
690 for(k = 0; k < order_e; ++k)
694 =
id * maxDof + orientMap[k] + 1;
711 m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
723 int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
727 int nTracePhys = trace->GetTotPoints();
737 nDim = locExpVector[0]->GetShapeDimension();
741 maxQuad = (1 > maxQuad ? 1 : maxQuad);
745 for (i = 0; i < trace->GetExpSize(); ++i)
747 quad = trace->GetExp(i)->GetTotPoints();
758 for (
int i = 0; i < trace->GetExpSize(); ++i)
760 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
761 offset = trace->GetPhys_Offset(i);
765 PeriodicMap::const_iterator it = perMap.find(eid);
766 if (perMap.count(eid) > 0)
771 eid = min(eid, ent.
id);
780 for (
int i = 0; i < trace->GetExpSize(); ++i)
782 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
783 offset = trace->GetPhys_Offset(i);
784 quad = trace->GetExp(i)->GetTotPoints();
790 PeriodicMap::const_iterator it = perMap.find(eid);
791 bool realign =
false;
792 if (perMap.count(eid) > 0)
797 realign = eid == min(eid, ent.
id);
798 eid = min(eid, ent.
id);
802 for (
int j = 0; j < quad; ++j)
813 it->second[0].orient, quad);
819 it->second[0].orient,
820 trace->GetExp(i)->GetNumPoints(0),
821 trace->GetExp(i)->GetNumPoints(1));
828 for (
int i = 0; i < nTracePhys; ++i)
834 for (
int i = 0; i < nTracePhys; ++i)
836 m_traceToUniversalMapUnique[i] = tmp2[i];
848 ASSERTL1(nquad2 == 0,
"nquad2 != 0 for reorienation");
849 for (
int i = 0; i < nquad1/2; ++i)
851 swap(toAlign[i], toAlign[nquad1-1-i]);
856 ASSERTL1(nquad2 != 0,
"nquad2 == 0 for reorienation");
866 for (
int i = 0; i < nquad2; ++i)
868 for (
int j = 0; j < nquad1; ++j)
870 tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
876 for (
int i = 0; i < nquad2; ++i)
878 for (
int j = 0; j < nquad1; ++j)
880 tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
891 for (
int i = 0; i < nquad2; ++i)
893 for (
int j = 0; j < nquad1/2; ++j)
895 swap(tmp[i*nquad1 + j],
896 tmp[i*nquad1 + nquad1-j-1]);
907 for (
int j = 0; j < nquad1; ++j)
909 for (
int i = 0; i < nquad2/2; ++i)
911 swap(tmp[i*nquad1 + j],
912 tmp[(nquad2-i-1)*nquad1 + j]);
1040 "i is out of range");
AssemblyMapDG()
Default constructor.
#define ASSERTL0(condition, msg)
bool m_signChange
Flag indicating if modes require sign reversal.
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
void MultiLevelBisectionReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm, BottomUpSubStructuredGraphSharedPtr &substructgraph, std::set< int > partVerts, int mdswitch)
LibUtilities::CommSharedPtr m_comm
Communicator.
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.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, int > m_bndCondTraceToGlobalTraceMap
Integer map of bnd cond trace number to global trace number.
bool isLocal
Flag specifying if this entity is local to this partition.
int GetBndSystemBandWidth() const
Returns the bandwidth of the boundary system.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Principle Modified Functions .
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
Lagrange Polynomials using the Gauss points .
Array< OneD, const NekDouble > GetLocalToGlobalBndSign() const
Retrieve the sign change for all local boundary modes.
int m_numLocalCoeffs
Total number of local coefficients.
void CuthillMckeeReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Array< OneD, int > m_traceToUniversalMap
Integer map of process trace space quadrature points to universal space.
std::vector< ExpansionSharedPtr > ExpansionVector
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Base class for constructing local to global mapping of degrees of freedom.
void UniversalTraceAssemble(Array< OneD, NekDouble > &pGlobal) const
size_t m_hash
Hash for map.
virtual void v_UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
const SpatialDomains::PointGeomSharedPtr GetGeom() const
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
boost::shared_ptr< SegExp > SegExpSharedPtr
Base class for all multi-elemental spectral/hp expansions.
Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > & GetElmtToTrace()
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique()
void SetUpUniversalTraceMap(const ExpList &locExp, const ExpListSharedPtr trace, const PeriodicMap &perMap=NullPeriodicMap)
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
GlobalSysSolnType m_solnType
The solution type of the global system.
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
int id
Geometry ID of entity.
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap()
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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.
boost::shared_ptr< PointExp > PointExpSharedPtr
void SetUpUniversalDGMap(const ExpList &locExp)
boost::shared_ptr< Expansion > ExpansionSharedPtr
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
int GetTraceToUniversalMap(int i)
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
int m_lowestStaticCondLevel
Lowest static condensation level.
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
virtual ~AssemblyMapDG()
Destructor.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Array< OneD, int > m_traceToUniversalMapUnique
Integer map of unique process trace space quadrature points to universal space (signed).
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
int GetTraceToUniversalMapUnique(int i)
int m_numLocalBndCoeffs
Number of local boundary coefficients.
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > m_elmtToTrace
list of edge expansions for a given element
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
virtual void v_GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
int GetOffset_Elmt_Id(int n) const
Get the element id associated with the n th consecutive block of data in m_phys and m_coeffs...
void AssembleBnd(const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign() const
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap()
SpatialDomains::Geometry1DSharedPtr GetGeom1D() const
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
int m_numGlobalCoeffs
Total number of global coefficients.
void GlobalToLocalBnd(const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
boost::shared_ptr< Geometry > GeometrySharedPtr
void RealignTraceElement(Array< OneD, int > &toAlign, StdRegions::Orientation orient, int nquad1, int nquad2=0)
Array< OneD, DataType > & GetPtr()
int GetNumDirichletBndPhys()
Return the number of boundary segments on which Dirichlet boundary conditions are imposed...
int m_numPatches
The number of patches (~elements) in the current level.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
int m_numDirichletBndPhys
Number of physical dirichlet boundary values in trace.
virtual int v_GetFullSystemBandWidth() const