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)
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)
588 locExpansion = locExpVector[eid];
589 nDim = locExpansion->GetShapeDimension();
594 int nverts = locExpansion->GetNverts();
595 for(j = 0; j < nverts; ++j)
599 id = locPointExp->
GetGeom()->GetGlobalID();
602 =
id * maxDof + j + 1;
608 for(j = 0; j < locExpansion->GetNedges(); ++j)
614 order_e = locExpansion->GetEdgeNcoeffs(j);
616 map<int,int> orientMap;
621 locExpansion->GetEdgeToElementMap(j, locExpansion->GetEorient(j), map2, sign2);
623 for (k = 0; k < map1.num_elements(); ++k)
628 for (
int l = 0; l < map2.num_elements(); ++l)
630 if (map1[k] == map2[l])
637 ASSERTL2(idx != -1,
"Problem with face to element map!");
641 for(k = 0; k < order_e; ++k)
645 =
id * maxDof + orientMap[k] + 1;
652 for(j = 0; j < locExpansion->GetNfaces(); ++j)
659 order_e = locExpansion->GetFaceNcoeffs(j);
661 map<int,int> orientMap;
666 locExpansion->GetFaceToElementMap(j, locExpansion->GetForient(j), map2, sign2);
668 for (k = 0; k < map1.num_elements(); ++k)
673 for (
int l = 0; l < map2.num_elements(); ++l)
675 if (map1[k] == map2[l])
682 ASSERTL2(idx != -1,
"Problem with face to element map!");
686 for(k = 0; k < order_e; ++k)
690 =
id * maxDof + orientMap[k] + 1;
707 m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
719 int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
723 int nTracePhys = trace->GetTotPoints();
733 nDim = locExpVector[0]->GetShapeDimension();
737 maxQuad = (1 > maxQuad ? 1 : maxQuad);
741 for (i = 0; i < trace->GetExpSize(); ++i)
743 quad = trace->GetExp(i)->GetTotPoints();
754 for (
int i = 0; i < trace->GetExpSize(); ++i)
756 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
757 offset = trace->GetPhys_Offset(i);
761 PeriodicMap::const_iterator it = perMap.find(eid);
762 if (perMap.count(eid) > 0)
767 eid = min(eid, ent.
id);
776 for (
int i = 0; i < trace->GetExpSize(); ++i)
778 eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
779 offset = trace->GetPhys_Offset(i);
780 quad = trace->GetExp(i)->GetTotPoints();
786 PeriodicMap::const_iterator it = perMap.find(eid);
787 bool realign =
false;
788 if (perMap.count(eid) > 0)
793 realign = eid == min(eid, ent.
id);
794 eid = min(eid, ent.
id);
798 for (
int j = 0; j < quad; ++j)
809 it->second[0].orient, quad);
815 it->second[0].orient,
816 trace->GetExp(i)->GetNumPoints(0),
817 trace->GetExp(i)->GetNumPoints(1));
824 for (
int i = 0; i < nTracePhys; ++i)
830 for (
int i = 0; i < nTracePhys; ++i)
832 m_traceToUniversalMapUnique[i] = tmp2[i];
844 ASSERTL1(nquad2 == 0,
"nquad2 != 0 for reorienation");
845 for (
int i = 0; i < nquad1/2; ++i)
847 swap(toAlign[i], toAlign[nquad1-1-i]);
852 ASSERTL1(nquad2 != 0,
"nquad2 == 0 for reorienation");
862 for (
int i = 0; i < nquad2; ++i)
864 for (
int j = 0; j < nquad1; ++j)
866 tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
872 for (
int i = 0; i < nquad2; ++i)
874 for (
int j = 0; j < nquad1; ++j)
876 tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
887 for (
int i = 0; i < nquad2; ++i)
889 for (
int j = 0; j < nquad1/2; ++j)
891 swap(tmp[i*nquad1 + j],
892 tmp[i*nquad1 + nquad1-j-1]);
903 for (
int j = 0; j < nquad1; ++j)
905 for (
int i = 0; i < nquad2/2; ++i)
907 swap(tmp[i*nquad1 + j],
908 tmp[(nquad2-i-1)*nquad1 + j]);
1038 "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.
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
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.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
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
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
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