53 CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(
60 const bool CheckforSingularSys):
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 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
253 ->GetGeom2D())->GetVid(j);
254 if(Dofs[0].count(vertId) == 0)
256 Dofs[0][vertId] = nvel*nz_loc;
259 if(IsDirVertDof.count(vertId) != 0)
261 Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc;
266 ->GetGeom2D())->GetEid(j);
267 if(Dofs[1].count(edgeId) == 0)
269 Dofs[1][edgeId] = nvel*(locExpVector[i]->GetEdgeNcoeffs(j)-2)*nz_loc;
273 if(IsDirEdgeDof.count(edgeId) != 0)
275 Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[i]->GetEdgeNcoeffs(j)-2);
280 set<int> extraDirVerts, extraDirEdges;
282 CreateGraph(*fields[0], bndCondExp, bndConditionsVec,
false,
283 periodicVerts, periodicEdges, periodicFaces,
284 ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
285 extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
310 edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
312 AddMeanPressureToEdgeId);
318 for(i = 0; i < nel; ++i)
320 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
323 ->GetGeom2D())->GetEid(j);
325 if(IsDirEdgeDof.count(edgeId) == 0)
329 if(AddMeanPressureToEdgeId[i] == -1)
331 AddMeanPressureToEdgeId[i] = edgeId;
335 ASSERTL0((AddMeanPressureToEdgeId[i] != -1),
"Did not determine "
336 "an edge to attach mean pressure dof");
338 Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
341 map<int,int> pressureEdgeOffset;
346 for(i = 0; i < bndCondExp.num_elements(); i++)
348 for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
350 bndSegExp = bndCondExp[i]->GetExp(j)
352 for(k = 0; k < nvel; ++k)
356 nLocDirBndCondDofs += bndSegExp->
GetNcoeffs()*nz_loc;
358 nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
383 Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
384 graphVertOffset[0] = 0;
388 for(i = 0; i < nel; ++i)
392 for(j = 0; j < locExpansion->GetNedges(); ++j)
396 ->GetGeom2D())->GetEid(j);
398 ->GetGeom2D())->GetVid(j);
400 for(k = 0; k < nvel*nz_loc; ++k)
402 graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+k] = 1;
403 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+k] = (nEdgeCoeffs-2);
406 bType = locExpansion->GetEdgeBasisType(j);
408 if( (nEdgeCoeffs >= 4)&&
418 for(i = 0; i < nel; ++i)
420 graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1] += nz_loc;
429 map<int,int> DirVertChk;
431 for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
434 for(j = 0; j < nvel; ++j)
443 if((cnt > 0)&&(cnt < nvel))
445 for(j = 0; j < nvel; ++j)
451 for(k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
454 id = bndCondExp[i]->GetExp(k)
456 ->GetGeom1D()->GetVid(0);
457 if(DirVertChk.count(
id*nvel+j) == 0)
459 DirVertChk[
id*nvel+j] = 1;
460 for(n = 0; n < nz_loc; ++n)
462 graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc + n] *= -1;
466 id = bndCondExp[i]->GetExp(k)
468 ->GetGeom1D()->GetVid(1);
469 if(DirVertChk.count(
id*nvel+j) == 0)
471 DirVertChk[
id*nvel+j] = 1;
472 for(n = 0; n < nz_loc; ++n)
474 graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc+n] *= -1;
479 id = bndCondExp[i]->GetExp(k)
481 ->GetGeom1D()->GetEid();
482 for(n = 0; n < nz_loc; ++n)
484 graphVertOffset[ReorderedGraphVertId[1][id]*nvel*nz_loc+j*nz_loc +n] *= -1;
495 for(i = 0; i < firstNonDirGraphVertId*nvel*nz_loc; ++i)
497 diff = abs(graphVertOffset[i]);
498 graphVertOffset[i] = cnt;
503 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
505 if(graphVertOffset[i] < 0)
507 diff = -graphVertOffset[i];
508 graphVertOffset[i] = -cnt;
517 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
519 if(graphVertOffset[i] >= 0)
521 diff = graphVertOffset[i];
522 graphVertOffset[i] = cnt;
529 for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
531 if(graphVertOffset[i] < 0)
533 graphVertOffset[i] = -graphVertOffset[i];
544 for(i = 0; i < nel; ++i)
571 for(i = 0; i < nel; ++i)
574 m_numLocalIntCoeffsPerPatch[i] = (
unsigned int) nz_loc*(pressure->GetExp(i)->GetNcoeffs()-1);
593 for(i = 0; i < nel; ++i)
603 map<int,int> inv_bmap;
604 locExpansion->GetBoundaryMap(bmap);
605 for(j = 0; j < bmap.num_elements(); ++j)
607 inv_bmap[bmap[j]] = j;
611 for(j = 0; j < locExpansion->GetNedges(); ++j)
613 nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
615 ->GetGeom2D())->GetEorient(j);
617 ->GetGeom2D())->GetEid(j);
619 ->GetGeom2D())->GetVid(j);
621 pIt = periodicEdges.find(meshEdgeId);
626 if (pIt != periodicEdges.end())
628 pair<int, StdRegions::Orientation> idOrient =
630 meshEdgeId, edgeOrient, pIt->second);
631 edgeOrient = idOrient.second;
634 locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
637 for(nv = 0; nv < nvel*nz_loc; ++nv)
639 m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
641 for(k = 0; k < nEdgeInteriorCoeffs; ++k)
643 m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] = graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
650 for(nv = 0; nv < nvel*nz_loc; ++nv)
652 for(k = 0; k < nEdgeInteriorCoeffs; ++k)
654 m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (
NekDouble) edgeInteriorSign[k];
661 nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc];
663 int psize = pressure->GetExp(i)->GetNcoeffs();
664 for(n = 0; n < nz_loc; ++n)
666 m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
668 pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
671 cnt += (velnbndry*nvel+ psize)*nz_loc;
676 for(nv = 0; nv < nvel; ++nv)
678 for(i = 0; i < bndCondExp.num_elements(); i++)
680 for(n = 0; n < nz_loc; ++n)
683 for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
685 bndSegExp = bndCondExp[i]->GetExp(j)
688 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
689 for(k = 0; k < 2; k++)
691 meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
692 m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->
GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
695 meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
697 nEdgeCoeffs = bndSegExp->GetNcoeffs();
698 for(k = 0; k < nEdgeCoeffs; k++)
700 if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
702 m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
703 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
707 ncoeffcnt += nEdgeCoeffs;
749 Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
750 for(i = 0; i < locExpVector.size(); ++i)
752 locExpansion = locExpVector[i]
754 for(j = 0; j < locExpansion->GetNverts(); ++j)
756 meshEdgeId = (locExpansion
758 ->GetGeom2D())->GetEid(j);
759 meshVertId = (locExpansion
761 ->GetGeom2D())->GetVid(j);
763 if(ReorderedGraphVertId[0][meshVertId] >=
764 firstNonDirGraphVertId)
766 vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
767 firstNonDirGraphVertId] =
771 if(ReorderedGraphVertId[1][meshEdgeId] >=
772 firstNonDirGraphVertId)
774 vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
775 firstNonDirGraphVertId] =
781 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
793 int &edgeId,
int &vertId,
int &firstNonDirGraphVertId, map<int,int> &IsDirEdgeDof,
802 map<int,int> HomGraphEdgeIdToEdgeId;
804 for(i = 0; i < nel; ++i)
806 for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
809 ->GetGeom2D())->GetEid(j);
812 if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
813 && (IsDirEdgeDof.count(edgeId) == 0))
815 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
817 if(EdgeIdToElmts[edgeId][0] == -1)
819 EdgeIdToElmts[edgeId][0] = i;
823 EdgeIdToElmts[edgeId][1] = i;
834 int nlevels = bottomUpGraph->GetNlevels();
840 vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
841 for(i = 0; i < bndgraphs.size(); ++i)
843 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
845 for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
848 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
850 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
865 for(
int n = 1; n < nlevels; ++n)
869 vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
872 map<int,int> ElmtInBndry;
874 for(i = 0; i < bndgraphs.size(); ++i)
876 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
878 for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
881 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
883 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
885 if(EdgeIdToElmts[edgeId][0] != -1)
887 ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
889 if(EdgeIdToElmts[edgeId][1] != -1)
891 ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
900 vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
901 for(i = 0; i < intgraphs.size(); ++i)
903 int GlobIdOffset = intgraphs[i]->GetIdOffset();
904 bool SetEdge =
false;
906 for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
909 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
911 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
913 for(k = 0; k < 2; ++k)
916 elmtid = EdgeIdToElmts[edgeId][k];
920 mapIt = ElmtInBndry.find(elmtid);
922 if(mapIt != ElmtInBndry.end())
925 int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
926 for(
int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
929 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
938 AddMeanPressureToEdgeId[elmtid] = defedge;
960 for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
962 if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
964 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
965 for(k = 0; k < 2; ++k)
968 elmtid = EdgeIdToElmts[edgeId][k];
981 if(AddMeanPressureToEdgeId[elmtid] == -1)
983 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)
Constructs mappings for the C0 scalar continuous Galerkin formulation.
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.