51 namespace MultiRegions
63 DisContField2D::DisContField2D(
void)
65 m_bndCondExpansions(),
73 const bool DeclareCoeffPhysArrays)
75 m_bndCondExpansions (In.m_bndCondExpansions),
76 m_bndConditions (In.m_bndConditions),
77 m_globalBndMat (In.m_globalBndMat),
78 m_traceMap (In.m_traceMap),
79 m_locTraceToTraceMap (In.m_locTraceToTraceMap),
80 m_boundaryEdges (In.m_boundaryEdges),
81 m_periodicVerts (In.m_periodicVerts),
82 m_periodicEdges (In.m_periodicEdges),
83 m_periodicFwdCopy (In.m_periodicFwdCopy),
84 m_periodicBwdCopy (In.m_periodicBwdCopy),
85 m_leftAdjacentEdges (In.m_leftAdjacentEdges)
90 *boost::dynamic_pointer_cast<ExpList1D>(In.
m_trace),
91 DeclareCoeffPhysArrays);
102 const std::string &variable,
103 const bool SetUpJustDG,
104 const bool DeclareCoeffPhysArrays)
105 :
ExpList2D(pSession, graph2D, DeclareCoeffPhysArrays, variable),
106 m_bndCondExpansions(),
115 if (variable.compare(
"DefaultVar") != 0)
119 DeclareCoeffPhysArrays);
121 if (DeclareCoeffPhysArrays)
147 for(e = 0; e < locExpList->GetExpSize(); ++e)
150 (*m_exp)[ElmtID[cnt+e]]->
151 as<LocalRegions::Expansion2D>();
153 locExpList->GetExp(e)->
154 as<LocalRegions::Expansion1D>();
156 locExpList->GetExp(e)->
157 as<LocalRegions::Expansion> ();
159 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
160 exp1d->SetAdjacentElementExp(EdgeID[cnt+e], exp2d);
165 if(
m_session->DefinesSolverInfo(
"PROJECTION"))
167 std::string ProjectStr =
169 if((ProjectStr ==
"MixedCGDG") ||
170 (ProjectStr ==
"Mixed_CG_Discontinuous"))
194 const std::string &variable,
195 const bool SetUpJustDG,
196 const bool DeclareCoeffPhysArrays)
202 if(variable.compare(
"DefaultVar") != 0)
207 if (DeclareCoeffPhysArrays)
235 for(e = 0; e < locExpList->GetExpSize(); ++e)
238 = (*m_exp)[ElmtID[cnt+e]]->
239 as<LocalRegions::Expansion2D>();
241 = locExpList->GetExp(e)->
242 as<LocalRegions::Expansion1D>();
244 = locExpList->GetExp(e)->
245 as<LocalRegions::Expansion> ();
247 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
248 exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
256 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
258 std::string ProjectStr =
261 if ((ProjectStr ==
"MixedCGDG") ||
262 (ProjectStr ==
"Mixed_CG_Discontinuous"))
317 for (e = 0; e < locExpList->GetExpSize(); ++e)
320 = (*m_exp)[ElmtID[cnt+e]]->
321 as<LocalRegions::Expansion2D>();
323 = locExpList->GetExp(e)->
324 as<LocalRegions::Expansion1D>();
326 = locExpList->GetExp(e)->
327 as<LocalRegions::Expansion> ();
329 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
330 exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
353 "Routine currently only tested for HybridDGHelmholtz");
356 "The local to global map is not set up for the requested "
365 (*m_globalBndMat)[mkey] = glo_matrix;
369 glo_matrix = matrixIter->second;
408 &elmtToTrace = m_traceMap->GetElmtToTrace();
415 for (
int i = 0; i <
m_exp->size(); ++i)
417 for (
int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
424 exp2d->SetEdgeExp (j, exp1d);
425 exp1d->SetAdjacentElementExp(j, exp2d);
433 for (
int i = 0; i <
m_trace->GetExpSize(); ++i)
438 int offset =
m_trace->GetPhys_Offset(i);
439 int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
444 if (traceGeomId != min(pIt->second[0].id, traceGeomId))
446 traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
447 traceEl->GetLeftAdjacentElementEdge());
450 else if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
452 traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
453 traceEl->GetLeftAdjacentElementEdge());
468 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
475 boost::unordered_map<int,pair<int,int> > perEdgeToExpMap;
476 boost::unordered_map<int,pair<int,int> >
::iterator it2;
477 for (cnt = n = 0; n <
m_exp->size(); ++n)
479 for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
482 (*
m_exp)[n]->GetGeom()->GetEid(e));
486 perEdgeToExpMap[it->first] = make_pair(n, e);
494 for (
int i = 0; i <
m_exp->size(); ++i)
496 for (
int j = 0; j < (*m_exp)[i]->GetNedges(); ++j, ++cnt)
504 for (
int n = 0; n <
m_exp->size(); ++n)
506 for (
int e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
508 int edgeGeomId = (*m_exp)[n]->GetGeom()->GetEid(e);
509 int offset =
m_trace->GetPhys_Offset(
510 elmtToTrace[n][e]->GetElmtId());
518 it2 = perEdgeToExpMap.find(ent.
id);
520 if (it2 == perEdgeToExpMap.end())
523 GetRowComm()->GetSize() > 1 && !ent.
isLocal)
529 ASSERTL1(
false,
"Periodic edge not found!");
534 "Periodic edge in non-forward space?");
536 int offset2 =
m_trace->GetPhys_Offset(
537 elmtToTrace[it2->second.first][it2->second.second]->
542 int nquad = elmtToTrace[n][e]->GetNumPoints(0);
544 vector<int> tmpBwd(nquad);
545 vector<int> tmpFwd(nquad);
549 for (
int i = 0; i < nquad; ++i)
551 tmpBwd[i] = offset2 + i;
552 tmpFwd[i] = offset + i;
557 for (
int i = 0; i < nquad; ++i)
559 tmpBwd[i] = offset2 + i;
560 tmpFwd[i] = offset + nquad - i - 1;
564 for (
int i = 0; i < nquad; ++i)
588 bool returnval =
true;
608 int vSame = (returnval?1:0);
632 const std::string &variable,
633 const bool DeclareCoeffPhysArrays)
642 SpatialDomains::BoundaryRegionCollection::const_iterator it;
645 for (it = bregions.begin(); it != bregions.end(); ++it)
663 for (it = bregions.begin(); it != bregions.end(); ++it)
671 DeclareCoeffPhysArrays, variable);
674 m_bndConditions[cnt] = bc;
677 std::string type = m_bndConditions[cnt]->GetUserDefined();
682 if((bc->GetBoundaryConditionType() !=
684 boost::iequals(type,
"I") ||
685 boost::iequals(type,
"CalcBC"))
708 const std::string &variable)
715 = boost::dynamic_pointer_cast<
717 SpatialDomains::BoundaryRegionCollection::const_iterator it;
731 map<int,int> perComps;
732 map<int,vector<int> > allVerts;
734 map<int,StdRegions::Orientation> allEdges;
736 int region1ID, region2ID, i, j, k, cnt;
740 for(i = 0; i < (*m_exp).size(); ++i)
742 for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
744 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
750 for (it = bregions.begin(); it != bregions.end(); ++it)
753 bconditions, it->first, variable);
755 if (locBCond->GetBoundaryConditionType()
762 region1ID = it->first;
763 region2ID = boost::static_pointer_cast<
765 locBCond)->m_connectedBoundaryRegion;
770 if (vComm->GetSize() == 1)
772 cId1 = it->second->begin()->first;
773 cId2 = bregions.find(region2ID)->second->begin()->first;
777 cId1 = bndRegOrder.find(region1ID)->second[0];
778 cId2 = bndRegOrder.find(region2ID)->second[0];
782 "Boundary region "+boost::lexical_cast<
string>(
783 region1ID)+
" should only contain 1 composite.");
788 vector<unsigned int> tmpOrder;
790 for (i = 0; i < c->size(); ++i)
793 boost::dynamic_pointer_cast<
795 ASSERTL0(segGeom,
"Unable to cast to shared ptr");
798 graph2D->GetElementsFromEdge(segGeom);
800 "The periodic boundaries belong to "
801 "more than one element of the mesh");
805 (*elmt)[0]->m_Element);
807 allEdges[(*c)[i]->GetGlobalID()] =
808 geom->GetEorient((*elmt)[0]->m_EdgeIndx);
812 if (vComm->GetSize() == 1)
814 tmpOrder.push_back((*c)[i]->GetGlobalID());
817 vector<int> vertList(2);
818 vertList[0] = segGeom->GetVid(0);
819 vertList[1] = segGeom->GetVid(1);
820 allVerts[(*c)[i]->GetGlobalID()] = vertList;
823 if (vComm->GetSize() == 1)
825 compOrder[it->second->begin()->first] = tmpOrder;
830 if (perComps.count(cId1) == 0)
832 if (perComps.count(cId2) == 0)
834 perComps[cId1] = cId2;
838 std::stringstream ss;
839 ss <<
"Boundary region " << cId2 <<
" should be "
840 <<
"periodic with " << perComps[cId2] <<
" but "
841 <<
"found " << cId1 <<
" instead!";
842 ASSERTL0(perComps[cId2] == cId1, ss.str());
847 std::stringstream ss;
848 ss <<
"Boundary region " << cId1 <<
" should be "
849 <<
"periodic with " << perComps[cId1] <<
" but "
850 <<
"found " << cId2 <<
" instead!";
851 ASSERTL0(perComps[cId1] == cId1, ss.str());
856 int n = vComm->GetSize();
857 int p = vComm->GetRank();
863 edgecounts[p] = allEdges.size();
867 for (i = 1; i < n; ++i)
869 edgeoffset[i] = edgeoffset[i-1] + edgecounts[i-1];
879 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
881 edgeIds [edgeoffset[p] + i ] = sIt->first;
882 edgeOrient[edgeoffset[p] + i ] = sIt->second;
883 edgeVerts [edgeoffset[p] + i++] = allVerts[sIt->first].size();
905 for (i = 0; i < n; ++i)
907 if (edgecounts[i] > 0)
910 edgecounts[i], edgeVerts + edgeoffset[i], 1);
919 for (i = 1; i < n; ++i)
921 vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
925 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
927 for (j = 0; j < allVerts[sIt->first].size(); ++j)
929 vertIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
936 map<int, StdRegions::Orientation> orientMap;
937 map<int, vector<int> > vertMap;
939 for (cnt = i = 0; i < totEdges; ++i)
941 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
942 "Already found edge in orientation map!");
945 vector<int> verts(edgeVerts[i]);
947 for (j = 0; j < edgeVerts[i]; ++j)
949 verts[j] = vertIds[cnt++];
951 vertMap[edgeIds[i]] = verts;
958 map<int,int>::const_iterator oIt;
960 map<int,int> allCompPairs;
969 for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
972 const int id1 = cIt->first;
973 const int id2 = cIt->second;
974 std::string id1s = boost::lexical_cast<
string>(id1);
975 std::string id2s = boost::lexical_cast<
string>(id2);
977 if (compMap.count(id1) > 0)
982 if (compMap.count(id2) > 0)
988 "Both composites not found on this process!");
993 map<int,int> compPairs;
996 "Unable to find composite "+id1s+
" in order map.");
998 "Unable to find composite "+id2s+
" in order map.");
999 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1000 "Periodic composites "+id1s+
" and "+id2s+
1001 " should have the same number of elements.");
1002 ASSERTL0(compOrder[id1].size() > 0,
1003 "Periodic composites "+id1s+
" and "+id2s+
1007 for (i = 0; i < compOrder[id1].size(); ++i)
1009 int eId1 = compOrder[id1][i];
1010 int eId2 = compOrder[id2][i];
1012 ASSERTL0(compPairs.count(eId1) == 0,
1015 if (compPairs.count(eId2) != 0)
1017 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1019 compPairs[eId1] = eId2;
1026 for (i = 0; i < 2; ++i)
1033 if (c[i]->size() > 0)
1035 for (j = 0; j < c[i]->size(); ++j)
1037 locEdges.insert(c[i]->at(j)->GetGlobalID());
1043 for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1045 int ids [2] = {pIt->first, pIt->second};
1046 bool local[2] = {locEdges.count(pIt->first) > 0,
1047 locEdges.count(pIt->second) > 0};
1049 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1050 orientMap.count(ids[1]) > 0,
1051 "Unable to find edge in orientation map.");
1053 allCompPairs[pIt->first ] = pIt->second;
1054 allCompPairs[pIt->second] = pIt->first;
1056 for (i = 0; i < 2; ++i)
1063 int other = (i+1) % 2;
1066 orientMap[ids[i]] == orientMap[ids[other]] ?
1075 for (i = 0; i < 2; ++i)
1077 int other = (i+1) % 2;
1080 orientMap[ids[i]] == orientMap[ids[other]] ?
1085 vector<int> perVerts1 = vertMap[ids[i]];
1086 vector<int> perVerts2 = vertMap[ids[other]];
1088 map<int, pair<int, bool> > tmpMap;
1092 tmpMap[perVerts1[0]] = make_pair(
1093 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1094 tmpMap[perVerts1[1]] = make_pair(
1095 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1099 tmpMap[perVerts1[0]] = make_pair(
1100 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1101 tmpMap[perVerts1[1]] = make_pair(
1102 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1105 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1110 mIt->second.second);
1114 if (perIt == periodicVerts.end())
1116 periodicVerts[mIt->first].push_back(ent2);
1117 perIt = periodicVerts.find(mIt->first);
1122 for (j = 0; j < perIt->second.size(); ++j)
1124 if (perIt->second[j].id == mIt->second.first)
1133 perIt->second.push_back(ent2);
1142 pairSizes[p] = allCompPairs.size();
1150 for (i = 1; i < n; ++i)
1152 pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1158 cnt = pairOffsets[p];
1160 for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1162 first [cnt ] = pIt->first;
1163 second[cnt++] = pIt->second;
1169 allCompPairs.clear();
1171 for(cnt = 0; cnt < totPairSizes; ++cnt)
1173 allCompPairs[first[cnt]] = second[cnt];
1179 for (cnt = i = 0; i < totEdges; ++i)
1181 int edgeId = edgeIds[i];
1183 ASSERTL0(allCompPairs.count(edgeId) > 0,
1184 "Unable to find matching periodic edge.");
1186 int perEdgeId = allCompPairs[edgeId];
1188 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1190 int vId = vertIds[cnt];
1194 if (perId == periodicVerts.end())
1201 orientMap[edgeId] == orientMap[perEdgeId] ?
1202 vertMap[perEdgeId][(j+1)%2] : vertMap[perEdgeId][j];
1206 locVerts.count(perVertexId) > 0);
1208 periodicVerts[vId].push_back(ent);
1216 for (perIt = periodicVerts.begin();
1217 perIt != periodicVerts.end(); ++perIt)
1220 for (i = 0; i < perIt->second.size(); ++i)
1222 perIt2 = periodicVerts.find(perIt->second[i].id);
1223 ASSERTL0(perIt2 != periodicVerts.end(),
1224 "Couldn't find periodic vertex.");
1226 for (j = 0; j < perIt2->second.size(); ++j)
1228 if (perIt2->second[j].id == perIt->first)
1234 for (k = 0; k < perIt->second.size(); ++k)
1236 if (perIt2->second[j].id == perIt->second[k].id)
1245 perIt->second.push_back(perIt2->second[j]);
1253 for (perIt = periodicVerts.begin();
1254 perIt != periodicVerts.end(); ++perIt)
1256 if (locVerts.count(perIt->first) > 0)
1268 as<LocalRegions::Expansion1D>();
1272 if (traceEl->GetLeftAdjacentElementEdge () == -1 ||
1273 traceEl->GetRightAdjacentElementEdge() == -1)
1283 int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
1289 fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1293 int offset =
m_trace->GetPhys_Offset(traceEl->GetElmtId());
1296 GetTraceToUniversalMapUnique(offset) >= 0;
1300 else if (traceEl->GetLeftAdjacentElementEdge () != -1 &&
1301 traceEl->GetRightAdjacentElementEdge() != -1)
1305 (traceEl->GetLeftAdjacentElementExp().get()) ==
1310 ASSERTL2(
false,
"Unconnected trace element!");
1357 int cnt, n, e,
npts, phys_offset;
1376 GetNFwdLocTracePts();
1385 boost::unordered_map<int,pair<int,int> >
::iterator it3;
1391 for(cnt = n = 0; n < nexp; ++n)
1396 for(e = 0; e < exp2d->GetNedges(); ++e, ++cnt)
1398 int offset =
m_trace->GetPhys_Offset(
1399 elmtToTrace[n][e]->GetElmtId());
1403 exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1404 field + phys_offset,
1405 e_tmp = Fwd + offset);
1409 exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1410 field + phys_offset,
1411 e_tmp = Bwd + offset);
1431 GetBndCondTraceToGlobalTraceMap(cnt+e));
1449 "Method not set up for non-zero Neumann "
1450 "boundary condition");
1451 id2 =
m_trace->GetPhys_Offset(
1452 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1462 "Method not set up for this boundary condition.");
1486 "Field must be in physical state to extract trace space.");
1508 int n, e, offset, phys_offset;
1514 "input array is of insufficient length");
1517 for(n = 0; n < nexp; ++n)
1521 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1523 offset =
m_trace->GetPhys_Offset(
1524 elmtToTrace[n][e]->GetElmtId());
1525 (*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
1526 inarray + phys_offset,
1527 e_tmp = outarray + offset);
1537 int e, n, offset, t_offset;
1545 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1547 t_offset =
GetTrace()->GetPhys_Offset(
1548 elmtToTrace[n][e]->GetElmtId());
1550 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1551 e,elmtToTrace[n][e],
1554 e_outarray = outarray+offset);
1585 m_trace->IProductWRTBase(Fn, Fcoeffs);
1592 int e, n, offset, t_offset;
1600 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1602 t_offset =
GetTrace()->GetPhys_Offset(
1603 elmtToTrace[n][e]->GetElmtId());
1604 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1605 e, elmtToTrace[n][e],
1607 e_outarray = outarray+offset);
1641 int e,n,offset, t_offset;
1649 for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1651 t_offset =
GetTrace()->GetPhys_Offset(
1652 elmtToTrace[n][e]->GetElmtId());
1657 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1658 e, elmtToTrace[n][e], Fwd+t_offset,
1659 e_outarray = outarray+offset);
1663 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1664 e, elmtToTrace[n][e], Bwd+t_offset,
1665 e_outarray = outarray+offset);
1680 map<int, int> globalIdMap;
1693 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1703 if(ElmtID.num_elements() != nbcs)
1708 if(EdgeID.num_elements() != nbcs)
1720 as<LocalRegions::Expansion1D>();
1723 graph2D->GetElementsFromEdge(exp1d->GetGeom1D());
1725 ElmtID[cnt] = globalIdMap[(*tmp)[0]->
1726 m_Element->GetGlobalID()];
1727 EdgeID[cnt] = (*tmp)[0]->m_EdgeIndx;
1733 boost::shared_ptr<ExpList> &result)
1736 int offsetOld, offsetNew;
1738 std::vector<unsigned int> eIDs;
1744 for (cnt = n = 0; n < i; ++n)
1752 eIDs.push_back(ElmtID[cnt+n]);
1760 for (n = 0; n < result->GetExpSize(); ++n)
1762 nq =
GetExp(ElmtID[cnt+n])->GetTotPoints();
1764 offsetNew = result->GetPhys_Offset(n);
1766 tmp2 = result->UpdatePhys()+ offsetNew, 1);
1768 nq =
GetExp(ElmtID[cnt+n])->GetNcoeffs();
1770 offsetNew = result->GetCoeff_Offset(n);
1772 tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
1804 int i,e,ncoeff_edge;
1814 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
1818 edge_lambda = loc_lambda;
1828 int nEdges = (*m_exp)[i]->GetNedges();
1831 for(e = 0; e < nEdges; ++e)
1833 edgedir = (*m_exp)[eid]->GetEorient(e);
1834 ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs();
1836 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
1837 elmtToTrace[eid][e]->SetCoeffsToOrientation(
1838 edgedir, edgeCoeffs[e], edgeCoeffs[e]);
1839 edge_lambda = edge_lambda + ncoeff_edge;
1842 (*m_exp)[eid]->DGDeriv(dir,
1846 out_tmp = out_d+cnt);
1847 cnt += (*m_exp)[eid]->GetNcoeffs();
1863 int i,j,n,cnt,cnt1,nbndry;
1880 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
1881 int NumDirichlet =
m_traceMap->GetNumLocalDirBndCoeffs();
1902 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
1911 for(cnt = cnt1 = n = 0; n < nexp; ++n)
1917 e_l = loc_lambda + cnt1;
1924 Floc =
Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1943 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1952 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1962 if(GloBndDofs - NumDirichlet > 0)
1981 m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1984 out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
2002 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2007 m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2008 LocLambda = (*HDGHelm) * LocLambda;
2009 m_traceMap->AssembleBnd(loc_lambda,outarray);
2027 map<int, RobinBCInfoSharedPtr> returnval;
2043 for(e = 0; e < locExpList->GetExpSize(); ++e)
2048 Array_tmp = locExpList->GetPhys() +
2049 locExpList->GetPhys_Offset(e));
2050 elmtid = ElmtID[cnt+e];
2052 if(returnval.count(elmtid) != 0)
2054 rInfo->next = returnval.find(elmtid)->second;
2056 returnval[elmtid] = rInfo;
2088 int i,cnt,e,ncoeff_edge;
2095 int eid,nq_elmt, nm_elmt;
2096 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2101 edge_lambda = loc_lambda;
2108 nq_elmt = (*m_exp)[eid]->GetTotPoints();
2109 nm_elmt = (*m_exp)[eid]->GetNcoeffs();
2113 out_tmp = force + nm_elmt;
2116 int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints();
2117 int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints();
2118 int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes();
2119 int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes();
2124 int nEdges = (*m_exp)[i]->GetNedges();
2127 for(e = 0; e < (*m_exp)[eid]->GetNedges(); ++e)
2129 edgedir = (*m_exp)[eid]->GetEorient(e);
2130 ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs();
2132 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
2133 elmtToTrace[eid][e]->SetCoeffsToOrientation(
2134 edgedir, edgeCoeffs[e], edgeCoeffs[e]);
2135 edge_lambda = edge_lambda + ncoeff_edge;
2163 ASSERTL0(
false,
"Wrong shape type, HDG postprocessing is not implemented");
2169 (*m_exp)[eid]->DGDeriv(
2171 elmtToTrace[eid], edgeCoeffs, out_tmp);
2172 (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2174 ppExp->IProductWRTDerivBase(0,qrhs,force);
2178 (*m_exp)[eid]->DGDeriv(
2180 elmtToTrace[eid], edgeCoeffs, out_tmp);
2182 (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2184 ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2189 (*m_exp)[eid]->BwdTrans(
2191 force[0] = (*m_exp)[eid]->Integral(qrhs);
2205 ppExp->BwdTrans(out.
GetPtr(), work);
2206 (*m_exp)[eid]->FwdTrans(work, tmp_coeffs = outarray +
m_coeff_offset[eid]);
2224 const std::string varName,
2234 for (i = 0; i < nbnd; ++i)
2240 npoints = locExpList->GetNpoints();
2248 locExpList->GetCoords(x0, x1, x2);
2252 locExpList->GetCoords(x0, x1, x2);
2259 string filebcs = boost::static_pointer_cast<
2270 boost::static_pointer_cast<
2273 m_dirichletCondition;
2275 condition.
Evaluate(x0, x1, x2, time,
2276 locExpList->UpdatePhys());
2279 locExpList->FwdTrans_BndConstrained(
2280 locExpList->GetPhys(),
2281 locExpList->UpdateCoeffs());
2286 string filebcs = boost::static_pointer_cast<
2297 boost::static_pointer_cast<
2301 condition.
Evaluate(x0, x1, x2, time,
2302 locExpList->UpdatePhys());
2305 locExpList->IProductWRTBase(
2306 locExpList->GetPhys(),
2307 locExpList->UpdateCoeffs());
2312 string filebcs = boost::static_pointer_cast<
2323 boost::static_pointer_cast<
2327 condition.
Evaluate(x0, x1, x2, time,
2328 locExpList->UpdatePhys());
2332 boost::static_pointer_cast<
2335 locExpList->IProductWRTBase(
2336 locExpList->GetPhys(),
2337 locExpList->UpdateCoeffs());
2342 locExpList->UpdatePhys());
2346 ASSERTL0(
false,
"This type of BC not implemented yet");
2356 locExpList->FwdTrans_IterPerExp(
2357 locExpList->GetPhys(),
2358 locExpList->UpdateCoeffs());
2362 ASSERTL0(
false,
"This type of BC not implemented yet");
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
virtual void v_FillBndCondFromField()
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
virtual void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This method extracts the trace (edges in 2D) from the field inarray and puts the values in outarray...
#define ASSERTL0(condition, msg)
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
void FindPeriodicEdges(const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
Determine the periodic edges and vertices for the given graph.
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
static ExpListSharedPtr NullExpListSharedPtr
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
bool isLocal
Flag specifying if this entity is local to this partition.
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
const BoundaryConditionCollection & GetBoundaryConditions(void) const
void ExtractFileBCs(const std::string &fileName, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
Lagrange Polynomials using the Gauss points .
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result)
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
std::map< int, std::vector< unsigned int > > BndRegionOrdering
std::map< ConstFactorType, NekDouble > ConstFactorMap
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
std::map< int, std::vector< unsigned int > > CompositeOrdering
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
boost::shared_ptr< MeshGraph2D > MeshGraph2DSharedPtr
int GetExpSize(void)
This function returns the number of elements in the expansion.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
Search through the edge expansions and identify which ones have Robin/Mixed type boundary conditions...
boost::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const Array< OneD, const NekDouble > &dirForcing)
Principle Orthogonal Functions .
bool SameTypeOfBoundaryConditions(const DisContField2D &In)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Base class for all multi-elemental spectral/hp expansions.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
DisContField2D()
Default constructor.
The base class for all shapes.
boost::shared_ptr< ExpList > & GetTrace()
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
std::vector< int > m_periodicBwdCopy
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
virtual ~DisContField2D()
Default destructor.
std::vector< bool > m_leftAdjacentEdges
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
boost::shared_ptr< SegGeom > SegGeomSharedPtr
std::set< int > m_boundaryEdges
A set storing the global IDs of any boundary edges.
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
NekDouble Evaluate() const
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Principle Orthogonal Functions .
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
LibUtilities::SessionReaderSharedPtr m_session
Session.
Defines a specification for a set of points.
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< Expansion > ExpansionSharedPtr
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
boost::shared_ptr< GeometryVector > Composite
bool IsLeftAdjacentEdge(const int n, const int e)
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
static const NekDouble kNekUnsetDouble
Describe a linear system.
std::map< int, Composite > CompositeMap
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
AssemblyMapDGSharedPtr m_traceMap
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::CommSharedPtr m_comm
Communicator.
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
GlobalLinSysMapShPtr m_globalBndMat
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
virtual void v_GetFwdBwdTracePhys(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
boost::shared_ptr< TriGeom > TriGeomSharedPtr
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
const BoundaryRegionCollection & GetBoundaryRegions(void) const
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Calculates the result of the multiplication of a global matrix of type specified by mkey with a vecto...
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
boost::shared_ptr< Basis > BasisSharedPtr
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph2D, const SpatialDomains::BoundaryConditions &bcs, const std::string &variable, const bool DeclareCoeffPhysArrays=true)
This function discretises the boundary conditions by setting up a list of one-dimensional boundary ex...
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Set up a list of element IDs and edge IDs that link to the boundary conditions.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error with respect to soln of the global spectral/hp element approximat...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
Describes the specification for a Basis.
Array< OneD, DataType > & GetPtr()
1D Gauss-Lobatto-Legendre quadrature points
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr