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>
58 namespace MultiRegions
61 m_numDirichletBndPhys(0)
77 const std::string variable):
80 int i, j, k, cnt, eid, id, id1, gid;
82 int nTraceExp = trace->GetExpSize();
83 int nbnd = bndCondExp.num_elements();
90 int nel = expList.size();
92 map<int, int> meshTraceId;
97 for(i = 0; i < nTraceExp; ++i)
99 meshTraceId[trace->GetExp(i)->GetGeom()->GetGlobalID()] = i;
104 for(i = 0; i < nel; ++i)
106 cnt += expList[i]->GetNtrace();
114 for(cnt = i = 0; i < nel; ++i)
118 for(j = 0; j < expList[i]->GetNtrace(); ++j)
120 id = expList[i]->GetGeom()->GetTid(j);
122 if(meshTraceId.count(
id) > 0)
125 trace->GetExp(meshTraceId.find(
id)->second);
129 ASSERTL0(
false,
"Failed to find trace map");
133 cnt += expList[i]->GetNtrace();
138 for(i = 0; i < nbnd; ++i)
140 cnt += bndCondExp[i]->GetExpSize();
149 for(i = 0; i < bndCondExp.num_elements(); ++i)
151 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
153 bndExp = bndCondExp[i]->GetExp(j);
154 traceGeom = bndExp->GetGeom();
155 id = traceGeom->GetGlobalID();
157 if(bndCond[i]->GetBoundaryConditionType() ==
178 for(i = 0; i < nel; ++i)
181 nbndry += expList[eid]->NumDGBndryCoeffs();
182 m_numLocalIntCoeffsPerPatch[i] = 0;
184 (
unsigned int) expList[eid]->NumDGBndryCoeffs();
200 typedef boost::adjacency_list<
201 boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
203 BoostGraph boostGraphObj;
204 int trace_id, trace_id1;
209 for(i = 0; i < nTraceExp; ++i)
211 id = trace->GetExp(i)->GetGeom()->GetGlobalID();
213 if (dirTrace.count(
id) == 0)
217 boost::add_vertex(boostGraphObj);
218 traceElmtGid[i] = cnt++;
223 traceElmtGid[i] = dirOffset;
224 dirOffset += trace->GetExp(i)->GetNcoeffs();
230 for(i = 0; i < nel; ++i)
234 for(j = 0; j < expList[eid]->GetNtrace(); ++j)
238 id = traceGeom->GetGlobalID();
239 trace_id = meshTraceId.find(
id)->second;
241 if(dirTrace.count(
id) == 0)
243 for(k = j+1; k < expList[eid]->GetNtrace(); ++k)
246 id1 = traceGeom->GetGlobalID();
247 trace_id1 = meshTraceId.find(id1)->second;
249 if(dirTrace.count(id1) == 0)
251 boost::add_edge((
size_t)traceElmtGid[trace_id],
252 (
size_t)traceElmtGid[trace_id1],
260 int nGraphVerts = nTraceExp - nDir;
266 for(i = 0; i < nGraphVerts; ++i)
268 vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
302 ASSERTL0(
false,
"Unrecognised solution type");
310 for(i = 0; i < nTraceExp - nDir; ++i)
312 traceElmtGid[perm[i]+nDir] = cnt;
313 cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
319 for(i = 0; i < nel; ++i)
325 for(j = 0; j < exp->GetNtrace(); ++j)
328 id = traceGeom->GetGlobalID();
329 gid = traceElmtGid[meshTraceId.find(
id)->second];
331 const int nDim = expList[eid]->GetNumBases();
340 order_e = expList[eid]->GetEdgeNcoeffs(j);
344 for(k = 0; k < order_e; ++k)
358 for (k = 2; k < order_e; ++k)
364 for(k = 3; k < order_e; k+=2)
366 m_localToGlobalBndSign[cnt+k] = -1.0;
373 for(k = 0; k < order_e; ++k)
382 for(k = 0; k < order_e; ++k)
390 ASSERTL0(
false,
"Boundary type not permitted");
397 order_e = expList[eid]->GetFaceNcoeffs(j);
399 std::map<int, int> orientMap;
410 expList[eid]->GetFaceToElementMap(j,fo,elmMap1,elmSign1);
411 expList[eid]->GetFaceToElementMap(
414 for (k = 0; k < elmMap1.num_elements(); ++k)
419 for (
int l = 0; l < elmMap2.num_elements(); ++l)
421 if (elmMap1[k] == elmMap2[l])
428 ASSERTL2(idx != -1,
"Problem with face to element map!");
432 for(k = 0; k < order_e; ++k)
435 m_localToGlobalBndSign[k+cnt] = elmSign2[orientMap[k]];
445 for(i = 0; i < nbnd; ++i)
447 cnt += bndCondExp[i]->GetNcoeffs();
453 int nbndexp = 0, bndOffset, bndTotal = 0;
454 for(cnt = i = 0; i < nbnd; ++i)
456 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
458 bndExp = bndCondExp[i]->GetExp(j);
459 id = bndExp->GetGeom()->GetGlobalID();
460 gid = traceElmtGid[meshTraceId.find(
id)->second];
461 bndOffset = bndCondExp[i]->GetCoeff_Offset(j) + bndTotal;
466 for(k = 0; k < bndExp->GetNcoeffs(); ++k)
472 nbndexp += bndCondExp[i]->GetExpSize();
473 bndTotal += bndCondExp[i]->GetNcoeffs();
491 for(
int i = 0; i < nGraphVerts; i++)
493 vwgts_perm[i] = vwgts[perm[i]];
496 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
504 for(i = 0; i < bndCondExp.num_elements(); ++i)
506 for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
508 bndExp = bndCondExp[i]->GetExp(j);
509 traceGeom = bndExp->GetGeom();
510 id = traceGeom->GetGlobalID();
512 meshTraceId.find(
id)->second;
554 for(i = 0; i < locExpVector.size(); ++i)
556 locExpansion = locExpVector[i];
557 nDim = locExpansion->GetShapeDimension();
562 maxDof = (1 > maxDof ? 1 : maxDof);
566 for (j = 0; j < locExpansion->GetNedges(); ++j)
568 dof = locExpansion->GetEdgeNcoeffs(j);
569 maxDof = (dof > maxDof ? dof : maxDof);
574 for (j = 0; j < locExpansion->GetNfaces(); ++j)
576 dof = locExpansion->GetFaceNcoeffs(j);
577 maxDof = (dof > maxDof ? dof : maxDof);
585 for(i = 0; i < locExpVector.size(); ++i)
587 locExpansion = locExpVector[i];
588 nDim = locExpansion->GetShapeDimension();
597 int nverts = locExpansion->GetNverts();
598 for(j = 0; j < nverts; ++j)
602 id = locPointExp->
GetGeom()->GetGlobalID();
605 =
id * maxDof + j + 1;
611 for(j = 0; j < locExpansion->GetNedges(); ++j)
617 order_e = locExpVector[eid]->GetEdgeNcoeffs(j);
619 map<int,int> orientMap;
624 locExpVector[eid]->GetEdgeToElementMap(j, locExpVector[eid]->GetEorient(j), map2, sign2);
626 for (k = 0; k < map1.num_elements(); ++k)
631 for (
int l = 0; l < map2.num_elements(); ++l)
633 if (map1[k] == map2[l])
640 ASSERTL2(idx != -1,
"Problem with face to element map!");
644 for(k = 0; k < order_e; ++k)
648 =
id * maxDof + orientMap[k] + 1;
655 for(j = 0; j < locExpansion->GetNfaces(); ++j)
662 order_e = locExpVector[eid]->GetFaceNcoeffs(j);
664 map<int,int> orientMap;
669 locExpVector[eid]->GetFaceToElementMap(j, locExpVector[eid]->GetForient(j), map2, sign2);
671 for (k = 0; k < map1.num_elements(); ++k)
676 for (
int l = 0; l < map2.num_elements(); ++l)
678 if (map1[k] == map2[l])
685 ASSERTL2(idx != -1,
"Problem with face to element map!");
689 for(k = 0; k < order_e; ++k)
693 =
id * maxDof + orientMap[k] + 1;
710 m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
722 int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
726 int nTracePhys = trace->GetTotPoints();
736 nDim = locExpVector[0]->GetShapeDimension();
740 maxQuad = (1 > maxQuad ? 1 : maxQuad);
744 for (i = 0; i < trace->GetExpSize(); ++i)
746 quad = trace->GetExp(i)->GetTotPoints();
757 for (
int i = 0; i < trace->GetExpSize(); ++i)
759 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
760 offset = trace->GetPhys_Offset(i);
764 PeriodicMap::const_iterator it = perMap.find(eid);
765 if (perMap.count(eid) > 0)
770 eid = min(eid, ent.
id);
779 for (
int i = 0; i < trace->GetExpSize(); ++i)
781 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
782 offset = trace->GetPhys_Offset(i);
783 quad = trace->GetExp(i)->GetTotPoints();
789 PeriodicMap::const_iterator it = perMap.find(eid);
790 bool realign =
false;
791 if (perMap.count(eid) > 0)
796 realign = eid == min(eid, ent.
id);
797 eid = min(eid, ent.
id);
801 for (
int j = 0; j < quad; ++j)
812 it->second[0].orient, quad);
818 it->second[0].orient,
819 trace->GetExp(i)->GetNumPoints(0),
820 trace->GetExp(i)->GetNumPoints(1));
827 for (
int i = 0; i < nTracePhys; ++i)
833 for (
int i = 0; i < nTracePhys; ++i)
835 m_traceToUniversalMapUnique[i] = tmp2[i];
847 ASSERTL1(nquad2 == 0,
"nquad2 != 0 for reorienation");
848 for (
int i = 0; i < nquad1/2; ++i)
850 swap(toAlign[i], toAlign[nquad1-1-i]);
855 ASSERTL1(nquad2 != 0,
"nquad2 == 0 for reorienation");
865 for (
int i = 0; i < nquad2; ++i)
867 for (
int j = 0; j < nquad1; ++j)
869 tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
875 for (
int i = 0; i < nquad2; ++i)
877 for (
int j = 0; j < nquad1; ++j)
879 tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
890 for (
int i = 0; i < nquad2; ++i)
892 for (
int j = 0; j < nquad1/2; ++j)
894 swap(tmp[i*nquad1 + j],
895 tmp[i*nquad1 + nquad1-j-1]);
906 for (
int j = 0; j < nquad1; ++j)
908 for (
int i = 0; i < nquad2/2; ++i)
910 swap(tmp[i*nquad1 + j],
911 tmp[(nquad2-i-1)*nquad1 + j]);
1039 "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.
std::map< int, vector< PeriodicEntity > > PeriodicMap
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)
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