58 const bool CheckforSingularSys)
62 size_t cnt = 0, offset = 0;
68 size_t nEdgeInteriorCoeffs;
69 int firstNonDirGraphVertId;
70 size_t nLocBndCondDofs = 0;
71 size_t nLocDirBndCondDofs = 0;
72 int nExtraDirichlet = 0;
79 size_t nvel = fields.size();
83 size_t nel = fields[0]->GetNumElmts();
88 std::vector<std::map<int, int>> ReorderedGraphVertId(3);
90 int staticCondLevel = 0;
92 if (CheckforSingularSys)
109 fields[0]->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
112 fields[0]->GetBndCondExpansions();
115 bndConditionsVec(nvel);
116 for (i = 0; i < nvel; ++i)
118 bndConditionsVec[i] = fields[i]->GetBndConditions();
121 std::map<int, int> IsDirVertDof;
122 std::map<int, int> IsDirEdgeDof;
125 for (j = 0; j < bndCondExp.size(); ++j)
127 std::map<int, int> BndExpVids;
129 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
135 BndExpVids[g->GetVid(0)] = g->GetVid(0);
136 BndExpVids[g->GetVid(1)] = g->GetVid(1);
139 for (i = 0; i < nvel; ++i)
141 if (bndConditionsVec[i][j]->GetBoundaryConditionType() ==
145 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
147 IsDirEdgeDof[bndCondExp[j]
151 ->GetGlobalID()] += 1;
157 for (
auto &mapIt : BndExpVids)
159 id = IsDirVertDof[mapIt.second] + 1;
160 IsDirVertDof[mapIt.second] = (
id > (int)nvel) ? nvel : id;
169 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
177 loc_exp->GetLeftAdjacentElementExp()->GetTraceNormal(
178 loc_exp->GetLeftAdjacentElementTrace());
180 size_t ndir = locnorm.size();
184 for (
size_t l = 0; l < locnorm[0].size(); ++l)
205 size_t edgeId, vertId;
214 for (i = 0; i < bndConditionsVec[0].size(); ++i)
216 if (bndConditionsVec[nvel - 1][i]->GetBoundaryConditionType() ==
228 ASSERTL0(
id != -1,
" Did not find an edge to attach singular pressure "
229 "degree of freedom");
233 for (i = 0; i < nel; ++i)
235 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
237 edgeId = (locExpVector[i]
242 if ((
int)edgeId == id)
244 AddMeanPressureToEdgeId[i] = id;
249 if (AddMeanPressureToEdgeId[i] != -1)
256 for (i = 0; i < nel; ++i)
258 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
263 if (Dofs[0].count(vertId) == 0)
265 Dofs[0][vertId] = nvel * nz_loc;
269 if (IsDirVertDof.count(vertId) != 0)
271 Dofs[0][vertId] -= IsDirVertDof[vertId] * nz_loc;
278 if (Dofs[1].count(edgeId) == 0)
281 nvel * (locExpVector[i]->GetTraceNcoeffs(j) - 2) * nz_loc;
286 if (IsDirEdgeDof.count(edgeId) != 0)
288 Dofs[1][edgeId] -= IsDirEdgeDof[edgeId] * nz_loc *
289 (locExpVector[i]->GetTraceNcoeffs(j) - 2);
294 std::set<int> extraDirVerts, extraDirEdges;
296 CreateGraph(*fields[0], bndCondExp, bndConditionsVec,
false, periodicVerts,
297 periodicEdges, periodicFaces, ReorderedGraphVertId,
298 bottomUpGraph, extraDirVerts, extraDirEdges,
299 firstNonDirGraphVertId, nExtraDirichlet, 4);
323 edgeId, vertId, firstNonDirGraphVertId,
324 IsDirEdgeDof, bottomUpGraph,
325 AddMeanPressureToEdgeId);
331 for (i = 0; i < nel; ++i)
333 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
339 if (IsDirEdgeDof.count(edgeId) == 0)
343 if (AddMeanPressureToEdgeId[i] == -1)
345 AddMeanPressureToEdgeId[i] = edgeId;
349 ASSERTL0((AddMeanPressureToEdgeId[i] != -1),
351 "an edge to attach mean pressure dof");
353 Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
356 std::map<int, int> pressureEdgeOffset;
361 for (i = 0; i < bndCondExp.size(); i++)
363 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
366 for (k = 0; k < nvel; ++k)
368 if (bndConditionsVec[k][i]->GetBoundaryConditionType() ==
371 nLocDirBndCondDofs += bndSegExp->
GetNcoeffs() * nz_loc;
374 if (bndConditionsVec[k][i]->GetBoundaryConditionType() !=
377 nLocBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
405 (ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),
407 graphVertOffset[0] = 0;
411 for (i = 0; i < nel; ++i)
415 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
418 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
419 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
421 for (k = 0; k < nvel * nz_loc; ++k)
423 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
426 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] * nvel *
428 k] = (nEdgeCoeffs - 2);
431 bType = locExpansion->GetBasisType(0);
441 for (i = 0; i < nel; ++i)
443 graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] +
456 std::map<int, int> DirVertChk;
458 for (i = 0; i < bndConditionsVec[0].size(); ++i)
461 for (j = 0; j < nvel; ++j)
463 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
471 if ((cnt > 0) && (cnt < nvel))
473 for (j = 0; j < nvel; ++j)
475 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
480 for (k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
488 if (DirVertChk.count(
id * nvel + j) == 0)
490 DirVertChk[
id * nvel + j] = 1;
491 for (n = 0; n < nz_loc; ++n)
493 graphVertOffset[ReorderedGraphVertId[0][id] *
495 j * nz_loc + n] *= -1;
504 if (DirVertChk.count(
id * nvel + j) == 0)
506 DirVertChk[
id * nvel + j] = 1;
507 for (n = 0; n < nz_loc; ++n)
509 graphVertOffset[ReorderedGraphVertId[0][id] *
511 j * nz_loc + n] *= -1;
521 for (n = 0; n < nz_loc; ++n)
523 graphVertOffset[ReorderedGraphVertId[1][id] * nvel *
525 j * nz_loc + n] *= -1;
535 for (i = 0; i < firstNonDirGraphVertId * nvel * nz_loc; ++i)
537 diff =
abs(graphVertOffset[i]);
538 graphVertOffset[i] = cnt;
543 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
546 if (graphVertOffset[i] < 0)
548 diff = -graphVertOffset[i];
549 graphVertOffset[i] = -cnt;
558 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
561 if (graphVertOffset[i] >= 0)
563 diff = graphVertOffset[i];
564 graphVertOffset[i] = cnt;
571 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
574 if (graphVertOffset[i] < 0)
576 graphVertOffset[i] = -graphVertOffset[i];
586 for (i = 0; i < nel; ++i)
589 nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
612 for (i = 0; i < nel; ++i)
615 (
unsigned int)nz_loc *
616 (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
618 (
unsigned int)nz_loc * (
pressure->GetExp(i)->GetNcoeffs() - 1);
629 for (i = 0; i < nel; ++i)
631 for (j = 0; j < nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs());
637 for (n = 0; n < nz_loc; ++n)
640 for (j = 1; j < (size_t)
pressure->GetExp(i)->GetNcoeffs(); ++j)
657 size_t nv, velnbndry;
662 for (i = 0; i < nel; ++i)
671 std::map<int, int> inv_bmap;
672 locExpansion->GetBoundaryMap(bmap);
673 for (j = 0; j < bmap.size(); ++j)
675 inv_bmap[bmap[j]] = j;
679 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
681 nEdgeInteriorCoeffs = locExpansion->GetTraceNcoeffs(j) - 2;
682 edgeOrient = (locExpansion->GetGeom())->GetEorient(j);
683 meshEdgeId = (locExpansion->GetGeom())->GetEid(j);
684 meshVertId = (locExpansion->GetGeom())->GetVid(j);
686 auto pIt = periodicEdges.find(meshEdgeId);
691 if (pIt != periodicEdges.end())
693 std::pair<int, StdRegions::Orientation> idOrient =
696 edgeOrient = idOrient.second;
699 locExpansion->GetTraceInteriorToElementMap(
700 j, edgeInteriorMap, edgeInteriorSign, edgeOrient);
703 for (nv = 0; nv < nvel * nz_loc; ++nv)
706 inv_bmap[locExpansion->GetVertexMap(j)]] =
707 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
712 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
715 inv_bmap[edgeInteriorMap[k]]] =
716 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] *
726 for (nv = 0; nv < nvel * nz_loc; ++nv)
728 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
731 inv_bmap[edgeInteriorMap[k]]] =
740 nEdgeInteriorCoeffs =
741 graphVertOffset[(ReorderedGraphVertId[1]
742 [AddMeanPressureToEdgeId[i]]) *
745 graphVertOffset[(ReorderedGraphVertId[1]
746 [AddMeanPressureToEdgeId[i]]) *
749 size_t psize =
pressure->GetExp(i)->GetNcoeffs();
750 for (n = 0; n < nz_loc; ++n)
754 [(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] + 1) *
757 nEdgeInteriorCoeffs +
758 pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
760 pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
763 cnt += (velnbndry * nvel + psize) * nz_loc;
768 for (nv = 0; nv < nvel; ++nv)
770 for (i = 0; i < bndCondExp.size(); i++)
772 if (bndConditionsVec[nv][i]->GetBoundaryConditionType() ==
778 for (n = 0; n < nz_loc; ++n)
781 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
786 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
787 for (k = 0; k < 2; k++)
789 meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
792 graphVertOffset[ReorderedGraphVertId[0]
798 meshEdgeId = (bndSegExp->GetGeom1D())->GetGlobalID();
800 nEdgeCoeffs = bndSegExp->GetNcoeffs();
801 for (k = 0; k < nEdgeCoeffs; k++)
807 [ReorderedGraphVertId[1][meshEdgeId] *
814 ncoeffcnt += nEdgeCoeffs;
857 firstNonDirGraphVertId);
858 for (i = 0; i < locExpVector.size(); ++i)
861 size_t nVert = locExpansion->
GetNverts();
862 for (j = 0; j < nVert; ++j)
864 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
865 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
867 if (ReorderedGraphVertId[0][meshVertId] >=
868 firstNonDirGraphVertId)
870 vwgts_perm[ReorderedGraphVertId[0][meshVertId] -
871 firstNonDirGraphVertId] =
875 if (ReorderedGraphVertId[1][meshEdgeId] >=
876 firstNonDirGraphVertId)
878 vwgts_perm[ReorderedGraphVertId[1][meshEdgeId] -
879 firstNonDirGraphVertId] =
885 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
895 std::vector<std::map<int, int>> &ReorderedGraphVertId,
size_t &nel,
897 [[maybe_unused]]
size_t &vertId,
int &firstNonDirGraphVertId,
898 std::map<int, int> &IsDirEdgeDof,
906 std::map<int, int> HomGraphEdgeIdToEdgeId;
908 for (i = 0; i < nel; ++i)
910 size_t nVert = locExpVector[i]->GetNverts();
911 for (j = 0; j < nVert; ++j)
918 if ((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId) &&
919 (IsDirEdgeDof.count(edgeId) == 0))
921 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId] -
922 firstNonDirGraphVertId] = edgeId;
924 if (EdgeIdToElmts[edgeId][0] == -1)
926 EdgeIdToElmts[edgeId][0] = i;
930 EdgeIdToElmts[edgeId][1] = i;
938 size_t nlevels = bottomUpGraph->GetNlevels();
944 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
945 bottomUpGraph->GetInteriorBlocks(nlevels);
946 for (i = 0; i < bndgraphs.size(); ++i)
948 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
949 size_t nVert = bndgraphs[i]->GetNverts();
951 for (j = 0; j < nVert; ++j)
954 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
956 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
971 for (
size_t n = 1; n < nlevels; ++n)
975 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
976 bottomUpGraph->GetInteriorBlocks(n + 1);
979 std::map<int, int> ElmtInBndry;
981 for (i = 0; i < bndgraphs.size(); ++i)
983 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
984 size_t nVert = bndgraphs[i]->GetNverts();
986 for (j = 0; j < nVert; ++j)
989 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
991 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
993 if (EdgeIdToElmts[edgeId][0] != -1)
995 ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
997 if (EdgeIdToElmts[edgeId][1] != -1)
999 ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
1008 std::vector<MultiRegions::SubGraphSharedPtr> intgraphs =
1009 bottomUpGraph->GetInteriorBlocks(n);
1010 for (i = 0; i < intgraphs.size(); ++i)
1012 int GlobIdOffset = intgraphs[i]->GetIdOffset();
1013 bool SetEdge =
false;
1015 size_t nVert = intgraphs[i]->GetNverts();
1016 for (j = 0; j < nVert; ++j)
1019 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1021 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1023 for (k = 0; k < 2; ++k)
1026 elmtid = EdgeIdToElmts[edgeId][k];
1030 auto mapIt = ElmtInBndry.find(elmtid);
1032 if (mapIt != ElmtInBndry.end())
1037 bndgraphs[mapIt->second]->GetIdOffset();
1039 l < bndgraphs[mapIt->second]->GetNverts();
1043 if (HomGraphEdgeIdToEdgeId.count(
1044 GlobIdOffset1 + l) != 0)
1046 AddMeanPressureToEdgeId[elmtid] =
1063 if (SetEdge ==
false)
1067 size_t nVert = intgraphs[i]->GetNverts();
1068 for (j = 0; j < nVert; ++j)
1070 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1072 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1073 for (k = 0; k < 2; ++k)
1076 elmtid = EdgeIdToElmts[edgeId][k];
1089 if (AddMeanPressureToEdgeId[elmtid] == -1)
1091 AddMeanPressureToEdgeId[elmtid] = defedge;
#define ASSERTL0(condition, msg)
CoupledLocalToGlobalC0ContMap(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const MultiRegions::ExpListSharedPtr &pressure, const size_t nz_loc, const bool CheeckForSingularSys=true)
void FindEdgeIdToAddMeanPressure(std::vector< std::map< int, int > > &ReorderedGraphVertId, size_t &nel, const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId, size_t &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Constructs mappings for the C0 scalar continuous Galerkin formulation.
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)
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numLocalCoeffs
Total number of local coefficients.
bool m_signChange
Flag indicating if modes require sign reversal.
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
int m_numGlobalCoeffs
Total number of global coefficients.
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
LibUtilities::SessionReaderSharedPtr m_session
Session object.
int m_numLocalBndCoeffs
Number of local boundary coefficients.
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
bool m_systemSingular
Flag indicating if the system is singular or not.
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
int m_numPatches
The number of patches (~elements) in the current level.
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Array< OneD, int > m_bndCondIDToGlobalTraceID
Integer map of bnd cond trace number to global trace number.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
int NumBndryCoeffs(void) const
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
int GetNverts() const
This function returns the number of vertices of the expansion domain.
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eModified_A
Principle Modified Functions .
std::shared_ptr< SegExp > SegExpSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
std::vector< ExpansionSharedPtr > ExpansionVector
@ eDirectMultiLevelStaticCond
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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...
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
static const NekDouble kNekZeroTol
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
scalarT< T > abs(scalarT< T > in)