53 CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(
60 const bool CheckforSingularSys):
61 AssemblyMapCG(pSession)
70 int nEdgeInteriorCoeffs;
71 int firstNonDirGraphVertId;
72 int nLocBndCondDofs = 0;
73 int nLocDirBndCondDofs = 0;
74 int nExtraDirichlet = 0;
81 int nvel = fields.num_elements();
82 MultiRegions::PeriodicMap::const_iterator pIt;
86 int nel = fields[0]->GetNumElmts();
91 vector<map<int,int> > ReorderedGraphVertId(3);
93 int staticCondLevel = 0;
95 if(CheckforSingularSys)
112 fields[0]->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
118 for(i = 0; i < nvel; ++i)
120 bndConditionsVec[i] = fields[i]->GetBndConditions();
123 map<int,int> IsDirVertDof;
124 map<int,int> IsDirEdgeDof;
128 for(j = 0; j < bndCondExp.num_elements(); ++j)
130 map<int,int> BndExpVids;
132 for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
136 BndExpVids[g->GetVid(0)] = g->GetVid(0);
137 BndExpVids[g->GetVid(1)] = g->GetVid(1);
140 for(i = 0; i < nvel; ++i)
145 for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
147 IsDirEdgeDof[bndCondExp[j]->GetExp(k)
149 ->GetGeom1D()->GetEid()] += 1;
156 for(mapIt = BndExpVids.begin(); mapIt != BndExpVids.end(); mapIt++)
158 id = IsDirVertDof[mapIt->second]+1;
159 IsDirVertDof[mapIt->second] = (
id > nvel)?nvel:
id;
168 for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
172 = bndCondExp[j]->GetExp(k)
174 locnorm = loc_exp->GetLeftAdjacentElementExp()->GetEdgeNormal(loc_exp->GetLeftAdjacentElementEdge());
177 int ndir = locnorm.num_elements();
180 for(
int l = 0; l < locnorm[0].num_elements(); ++l)
211 for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
215 id = bndCondExp[i]->GetExp(0)
222 ASSERTL0(
id != -1,
" Did not find an edge to attach singular pressure degree of freedom");
226 for(i = 0; i < nel; ++i)
228 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
231 ->GetGeom2D())->GetEid(j);
235 AddMeanPressureToEdgeId[i] = id;
240 if(AddMeanPressureToEdgeId[i] != -1)
248 for(i = 0; i < nel; ++i)
250 eid = fields[0]->GetOffset_Elmt_Id(i);
251 for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
254 ->GetGeom2D())->GetVid(j);
255 if(Dofs[0].count(vertId) == 0)
257 Dofs[0][vertId] = nvel*nz_loc;
260 if(IsDirVertDof.count(vertId) != 0)
262 Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc;
267 ->GetGeom2D())->GetEid(j);
268 if(Dofs[1].count(edgeId) == 0)
270 Dofs[1][edgeId] = nvel*(locExpVector[eid]->GetEdgeNcoeffs(j)-2)*nz_loc;
274 if(IsDirEdgeDof.count(edgeId) != 0)
276 Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[eid]->GetEdgeNcoeffs(j)-2);
281 set<int> extraDirVerts, extraDirEdges;
283 CreateGraph(*fields[0], bndCondExp, bndConditionsVec,
false,
284 periodicVerts, periodicEdges, periodicFaces,
285 ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
286 extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
311 edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
313 AddMeanPressureToEdgeId);
319 for(i = 0; i < nel; ++i)
321 eid = fields[0]->GetOffset_Elmt_Id(i);
322 for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
325 ->GetGeom2D())->GetEid(j);
327 if(IsDirEdgeDof.count(edgeId) == 0)
331 if(AddMeanPressureToEdgeId[eid] == -1)
333 AddMeanPressureToEdgeId[eid] = edgeId;
337 ASSERTL0((AddMeanPressureToEdgeId[eid] != -1),
"Did not determine "
338 "an edge to attach mean pressure dof");
340 Dofs[1][AddMeanPressureToEdgeId[eid]] += nz_loc;
343 map<int,int> pressureEdgeOffset;
348 for(i = 0; i < bndCondExp.num_elements(); i++)
350 for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
352 bndSegExp = bndCondExp[i]->GetExp(j)
354 for(k = 0; k < nvel; ++k)
358 nLocDirBndCondDofs += bndSegExp->
GetNcoeffs()*nz_loc;
360 nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
385 Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
386 graphVertOffset[0] = 0;
390 for(i = 0; i < nel; ++i)
392 eid = fields[0]->GetOffset_Elmt_Id(i);
395 for(j = 0; j < locExpansion->GetNedges(); ++j)
399 ->GetGeom2D())->GetEid(j);
401 ->GetGeom2D())->GetVid(j);
403 for(k = 0; k < nvel*nz_loc; ++k)
405 graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+k] = 1;
406 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+k] = (nEdgeCoeffs-2);
409 bType = locExpansion->GetEdgeBasisType(j);
411 if( (nEdgeCoeffs >= 4)&&
421 for(i = 0; i < nel; ++i)
423 graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1] += nz_loc;
432 map<int,int> DirVertChk;
434 for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
437 for(j = 0; j < nvel; ++j)
446 if((cnt > 0)&&(cnt < nvel))
448 for(j = 0; j < nvel; ++j)
454 for(k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
457 id = bndCondExp[i]->GetExp(k)
459 ->GetGeom1D()->GetVid(0);
460 if(DirVertChk.count(
id*nvel+j) == 0)
462 DirVertChk[
id*nvel+j] = 1;
463 for(n = 0; n < nz_loc; ++n)
465 graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc + n] *= -1;
469 id = bndCondExp[i]->GetExp(k)
471 ->GetGeom1D()->GetVid(1);
472 if(DirVertChk.count(
id*nvel+j) == 0)
474 DirVertChk[
id*nvel+j] = 1;
475 for(n = 0; n < nz_loc; ++n)
477 graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc+n] *= -1;
482 id = bndCondExp[i]->GetExp(k)
484 ->GetGeom1D()->GetEid();
485 for(n = 0; n < nz_loc; ++n)
487 graphVertOffset[ReorderedGraphVertId[1][id]*nvel*nz_loc+j*nz_loc +n] *= -1;
498 for(i = 0; i < firstNonDirGraphVertId*nvel*nz_loc; ++i)
500 diff = abs(graphVertOffset[i]);
501 graphVertOffset[i] = cnt;
506 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
508 if(graphVertOffset[i] < 0)
510 diff = -graphVertOffset[i];
511 graphVertOffset[i] = -cnt;
520 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
522 if(graphVertOffset[i] >= 0)
524 diff = graphVertOffset[i];
525 graphVertOffset[i] = cnt;
532 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
534 if(graphVertOffset[i] < 0)
536 graphVertOffset[i] = -graphVertOffset[i];
547 for(i = 0; i < nel; ++i)
574 for(i = 0; i < nel; ++i)
577 m_numLocalIntCoeffsPerPatch[i] = (
unsigned int) nz_loc*(pressure->GetExp(i)->GetNcoeffs()-1);
596 for(i = 0; i < nel; ++i)
598 eid = fields[0]->GetOffset_Elmt_Id(i);
607 map<int,int> inv_bmap;
608 locExpansion->GetBoundaryMap(bmap);
609 for(j = 0; j < bmap.num_elements(); ++j)
611 inv_bmap[bmap[j]] = j;
615 for(j = 0; j < locExpansion->GetNedges(); ++j)
617 nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
619 ->GetGeom2D())->GetEorient(j);
621 ->GetGeom2D())->GetEid(j);
623 ->GetGeom2D())->GetVid(j);
625 pIt = periodicEdges.find(meshEdgeId);
630 if (pIt != periodicEdges.end())
632 pair<int, StdRegions::Orientation> idOrient =
634 meshEdgeId, edgeOrient, pIt->second);
635 edgeOrient = idOrient.second;
638 locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
641 for(nv = 0; nv < nvel*nz_loc; ++nv)
643 m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
645 for(k = 0; k < nEdgeInteriorCoeffs; ++k)
647 m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] = graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
654 for(nv = 0; nv < nvel*nz_loc; ++nv)
656 for(k = 0; k < nEdgeInteriorCoeffs; ++k)
658 m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (
NekDouble) edgeInteriorSign[k];
665 nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc];
667 int psize = pressure->GetExp(eid)->GetNcoeffs();
668 for(n = 0; n < nz_loc; ++n)
670 m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[eid]];
672 pressureEdgeOffset[AddMeanPressureToEdgeId[eid]] += 1;
675 cnt += (velnbndry*nvel+ psize)*nz_loc;
680 for(nv = 0; nv < nvel; ++nv)
682 for(i = 0; i < bndCondExp.num_elements(); i++)
684 for(n = 0; n < nz_loc; ++n)
687 for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
689 bndSegExp = bndCondExp[i]->GetExp(j)
692 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
693 for(k = 0; k < 2; k++)
695 meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
696 m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->
GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
699 meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
701 nEdgeCoeffs = bndSegExp->GetNcoeffs();
702 for(k = 0; k < nEdgeCoeffs; k++)
704 if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
706 m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
707 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
711 ncoeffcnt += nEdgeCoeffs;
753 Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
754 for(i = 0; i < locExpVector.size(); ++i)
756 eid = fields[0]->GetOffset_Elmt_Id(i);
757 locExpansion = locExpVector[eid]
759 for(j = 0; j < locExpansion->GetNverts(); ++j)
761 meshEdgeId = (locExpansion
763 ->GetGeom2D())->GetEid(j);
764 meshVertId = (locExpansion
766 ->GetGeom2D())->GetVid(j);
768 if(ReorderedGraphVertId[0][meshVertId] >=
769 firstNonDirGraphVertId)
771 vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
772 firstNonDirGraphVertId] =
776 if(ReorderedGraphVertId[1][meshEdgeId] >=
777 firstNonDirGraphVertId)
779 vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
780 firstNonDirGraphVertId] =
786 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
798 int &edgeId,
int &vertId,
int &firstNonDirGraphVertId, map<int,int> &IsDirEdgeDof,
807 map<int,int> HomGraphEdgeIdToEdgeId;
809 for(i = 0; i < nel; ++i)
811 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
814 ->GetGeom2D())->GetEid(j);
817 if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
818 && (IsDirEdgeDof.count(edgeId) == 0))
820 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
822 if(EdgeIdToElmts[edgeId][0] == -1)
824 EdgeIdToElmts[edgeId][0] = i;
828 EdgeIdToElmts[edgeId][1] = i;
839 int nlevels = bottomUpGraph->GetNlevels();
845 vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
846 for(i = 0; i < bndgraphs.size(); ++i)
848 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
850 for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
853 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
855 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
870 for(
int n = 1; n < nlevels; ++n)
874 vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
877 map<int,int> ElmtInBndry;
879 for(i = 0; i < bndgraphs.size(); ++i)
881 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
883 for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
886 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
888 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
890 if(EdgeIdToElmts[edgeId][0] != -1)
892 ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
894 if(EdgeIdToElmts[edgeId][1] != -1)
896 ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
905 vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
906 for(i = 0; i < intgraphs.size(); ++i)
908 int GlobIdOffset = intgraphs[i]->GetIdOffset();
909 bool SetEdge =
false;
911 for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
914 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
916 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
918 for(k = 0; k < 2; ++k)
921 elmtid = EdgeIdToElmts[edgeId][k];
925 mapIt = ElmtInBndry.find(elmtid);
927 if(mapIt != ElmtInBndry.end())
930 int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
931 for(
int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
934 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
943 AddMeanPressureToEdgeId[elmtid] = defedge;
965 for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
967 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
969 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
970 for(k = 0; k < 2; ++k)
973 elmtid = EdgeIdToElmts[edgeId][k];
986 if(AddMeanPressureToEdgeId[elmtid] == -1)
988 AddMeanPressureToEdgeId[elmtid] = defedge;
#define ASSERTL0(condition, msg)
bool m_systemSingular
Flag indicating if the system is singular or not.
bool m_signChange
Flag indicating if modes require sign reversal.
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
Principle Modified Functions .
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
int m_numLocalCoeffs
Total number of local coefficients.
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< ExpansionSharedPtr > ExpansionVector
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
int NumBndryCoeffs(void) const
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
static const NekDouble kNekZeroTol
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Principle Modified Functions .
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
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, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
LibUtilities::SessionReaderSharedPtr m_session
Session object.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
void FindEdgeIdToAddMeanPressure(std::vector< std::map< int, int > > &ReorderedGraphVertId, int &nel, const LocalRegions::ExpansionVector &locExpVector, int &edgeId, int &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
int m_numGlobalCoeffs
Total number of global coefficients.
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
int m_numPatches
The number of patches (~elements) in the current level.