63 :
ExpList(), m_bndConditions(), m_bndCondExpansions(),
82 const std::string &variable,
const bool SetUpJustDG,
83 const bool DeclareCoeffPhysArrays,
85 const std::string bcvariable)
86 :
ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
90 if (bcvariable ==
"NotSet")
99 if (bcvar.compare(
"DefaultVar") != 0)
104 DeclareCoeffPhysArrays);
105 if (DeclareCoeffPhysArrays)
122 for (
int e = 0; e < locExpList->GetExpSize(); ++e)
124 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
125 locExpList->GetExp(e));
126 locExpList->GetExp(e)->SetAdjacentElementExp(
127 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
138 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
140 std::string ProjectStr =
m_session->GetSolverInfo(
"PROJECTION");
141 if ((ProjectStr ==
"MixedCGDG") ||
142 (ProjectStr ==
"Mixed_CG_Discontinuous"))
176 m_comm,
true,
"DefaultVar", ImpType);
186 if (
m_session->DefinesCmdLineArgument(
"verbose"))
199 for (
int i = 0; i <
m_exp->size(); ++i)
201 for (
int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
206 exp->SetAdjacentElementExp(j, (*
m_exp)[i]);
207 (*m_exp)[i]->SetTraceExp(j, exp);
229 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
236 std::unordered_map<int, pair<int, int>> perTraceToExpMap;
237 for (cnt = n = 0; n <
m_exp->size(); ++n)
239 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
241 auto it = periodicTraces.find((*
m_exp)[n]->GetGeom()->GetTid(v));
243 if (it != periodicTraces.end())
245 perTraceToExpMap[it->first] = make_pair(n, v);
254 for (cnt = n = 0; n <
m_exp->size(); ++n)
256 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
263 for (cnt = n = 0; n <
m_exp->size(); ++n)
265 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
267 int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
270 auto it = periodicTraces.find(GeomId);
272 if (it != periodicTraces.end())
275 auto it2 = perTraceToExpMap.find(ent.
id);
277 if (it2 == perTraceToExpMap.end())
279 if (
m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
286 ASSERTL1(
false,
"Periodic trace not found!");
291 "Periodic trace in non-forward space?");
293 int offset =
m_trace->GetPhys_Offset(
294 (
m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
296 int offset2 =
m_trace->GetPhys_Offset(
297 (
m_traceMap->GetElmtToTrace())[it2->second.first]
313 int nquad = elmtToTrace[n][v]->GetNumPoints(0);
315 vector<int> tmpBwd(nquad);
316 vector<int> tmpFwd(nquad);
320 for (
int i = 0; i < nquad; ++i)
322 tmpBwd[i] = offset2 + i;
323 tmpFwd[i] = offset + i;
328 for (
int i = 0; i < nquad; ++i)
330 tmpBwd[i] = offset2 + i;
331 tmpFwd[i] = offset + nquad - i - 1;
335 for (
int i = 0; i < nquad; ++i)
346 int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
347 int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
349 vector<int> tmpBwd(nquad1 * nquad2);
350 vector<int> tmpFwd(nquad1 * nquad2);
360 for (
int i = 0; i < nquad2; ++i)
362 for (
int j = 0; j < nquad1; ++j)
364 tmpBwd[i * nquad1 + j] =
365 offset2 + i * nquad1 + j;
366 tmpFwd[i * nquad1 + j] =
367 offset + j * nquad2 + i;
373 for (
int i = 0; i < nquad2; ++i)
375 for (
int j = 0; j < nquad1; ++j)
377 tmpBwd[i * nquad1 + j] =
378 offset2 + i * nquad1 + j;
379 tmpFwd[i * nquad1 + j] =
380 offset + i * nquad1 + j;
394 for (
int i = 0; i < nquad2; ++i)
396 for (
int j = 0; j < nquad1 / 2; ++j)
398 swap(tmpFwd[i * nquad1 + j],
399 tmpFwd[i * nquad1 + nquad1 - j - 1]);
413 for (
int j = 0; j < nquad1; ++j)
415 for (
int i = 0; i < nquad2 / 2; ++i)
417 swap(tmpFwd[i * nquad1 + j],
418 tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
423 for (
int i = 0; i < nquad1 * nquad2; ++i)
486 if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
487 traceEl->GetRightAdjacentElementTrace() == -1)
500 else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
501 traceEl->GetRightAdjacentElementTrace() != -1)
504 fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*
m_exp)[n].get();
508 ASSERTL2(
false,
"Unconnected trace element!");
520 [[maybe_unused]]
const std::string &variable)
527 map<int, int> GeometryToRegionsMap;
535 for (
auto &it : bregions)
537 for (
auto &bregionIt : *it.second)
541 int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
542 GeometryToRegionsMap[id] = it.first;
546 map<int, SpatialDomains::GeometrySharedPtr> EndOfDomain;
549 for (
auto &domIt : domain)
552 for (
int i = 0; i < geomvector->m_geomVec.size(); ++i)
554 for (
int j = 0; j < 2; ++j)
556 int vid = geomvector->m_geomVec[i]->GetVid(j);
557 if (EndOfDomain.count(vid) == 0)
559 EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
563 EndOfDomain.erase(vid);
568 ASSERTL1(EndOfDomain.size() == 2,
"Did not find two ends of domain");
570 for (
auto ®It : EndOfDomain)
572 if (GeometryToRegionsMap.count(regIt.first) != 0)
575 auto iter = GeometryToRegionsMap.find(regIt.first);
576 ASSERTL1(iter != GeometryToRegionsMap.end(),
577 "Failied to find GeometryToRegionMap");
579 int regionId = iter->second;
580 auto bregionsIter = bregions.find(regionId);
581 ASSERTL1(bregionsIter != bregions.end(),
582 "Failed to find boundary region");
585 returnval->AddBoundaryRegions(regionId, breg);
587 auto bconditionsIter = bconditions.find(regionId);
588 ASSERTL1(bconditionsIter != bconditions.end(),
589 "Failed to find boundary collection");
592 bconditionsIter->second;
593 returnval->AddBoundaryConditions(regionId, bcond);
611 const std::string &variable,
613 bool SetToOneSpaceDimension,
615 :
ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
618 if (variable.compare(
"DefaultVar") != 0)
637 const bool DeclareCoeffPhysArrays)
638 :
ExpList(In, DeclareCoeffPhysArrays), m_bndConditions(In.m_bndConditions),
639 m_bndCondExpansions(In.m_bndCondExpansions),
640 m_globalBndMat(In.m_globalBndMat), m_traceMap(In.m_traceMap),
641 m_boundaryTraces(In.m_boundaryTraces),
642 m_periodicVerts(In.m_periodicVerts),
643 m_periodicFwdCopy(In.m_periodicFwdCopy),
644 m_periodicBwdCopy(In.m_periodicBwdCopy),
645 m_leftAdjacentTraces(In.m_leftAdjacentTraces),
646 m_locTraceToTraceMap(In.m_locTraceToTraceMap)
651 *In.
m_trace, DeclareCoeffPhysArrays);
661 const std::string &variable,
const bool SetUpJustDG,
662 const bool DeclareCoeffPhysArrays)
663 :
ExpList(In, DeclareCoeffPhysArrays)
670 if (variable.compare(
"DefaultVar") != 0)
675 if (DeclareCoeffPhysArrays)
703 for (e = 0; e < locExpList->GetExpSize(); ++e)
705 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
706 TraceID[cnt + e], locExpList->GetExp(e));
707 locExpList->GetExp(e)->SetAdjacentElementExp(
708 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
714 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
716 std::string ProjectStr =
719 if ((ProjectStr ==
"MixedCGDG") ||
720 (ProjectStr ==
"Mixed_CG_Discontinuous"))
764 for (e = 0; e < locExpList->GetExpSize(); ++e)
766 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
767 TraceID[cnt + e], locExpList->GetExp(e));
768 locExpList->GetExp(e)->SetAdjacentElementExp(
769 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
815 const bool DeclareCoeffPhysArrays)
833 for (
auto &it : bregions)
838 m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
839 false, bc->GetComm());
850 boost::iequals(type,
"I") || boost::iequals(type,
"CalcBC"))
879 int i, region1ID, region2ID;
882 map<int, int> BregionToVertMap;
886 for (
auto &it : bregions)
891 if (locBCond->GetBoundaryConditionType() !=
897 it.second->begin()->second->m_geomVec[0]->GetGlobalID();
899 BregionToVertMap[it.first] = id;
904 int n = vComm->GetSize();
905 int p = vComm->GetRank();
908 nregions[
p] = BregionToVertMap.size();
915 for (i = 1; i < n; ++i)
917 regOffset[i] = regOffset[i - 1] + nregions[i - 1];
924 for (
auto &iit : BregionToVertMap)
926 bregid[i] = iit.first;
927 bregmap[i++] = iit.second;
928 islocal.insert(iit.first);
934 for (
int i = 0; i < totRegions; ++i)
936 BregionToVertMap[bregid[i]] = bregmap[i];
940 for (
auto &it : bregions)
945 if (locBCond->GetBoundaryConditionType() !=
952 region1ID = it.first;
954 std::static_pointer_cast<
956 ->m_connectedBoundaryRegion;
958 ASSERTL0(BregionToVertMap.count(region1ID) != 0,
959 "Cannot determine vertex of region1ID");
961 ASSERTL0(BregionToVertMap.count(region2ID) != 0,
962 "Cannot determine vertex of region2ID");
966 islocal.count(region2ID) != 0);
974 int region1ID, region2ID, i, j, k, cnt;
978 m_graph->GetCompositeOrdering();
980 m_graph->GetBndRegionOrdering();
987 map<int, int> perComps;
988 map<int, vector<int>> allVerts;
990 map<int, pair<int, StdRegions::Orientation>> allEdges;
993 for (i = 0; i < (*m_exp).size(); ++i)
995 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
997 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1003 for (
auto &it : bregions)
1008 if (locBCond->GetBoundaryConditionType() !=
1015 region1ID = it.first;
1017 std::static_pointer_cast<
1019 ->m_connectedBoundaryRegion;
1024 if (vComm->GetSize() == 1)
1026 cId1 = it.second->begin()->first;
1027 cId2 = bregions.find(region2ID)->second->begin()->first;
1031 cId1 = bndRegOrder.find(region1ID)->second[0];
1032 cId2 = bndRegOrder.find(region2ID)->second[0];
1036 "Boundary region " + std::to_string(region1ID) +
1037 " should only contain 1 composite.");
1039 vector<unsigned int> tmpOrder;
1044 it.second->begin()->second;
1046 for (i = 0; i < c->m_geomVec.size(); ++i)
1049 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1051 ASSERTL0(segGeom,
"Unable to cast to shared ptr");
1054 m_graph->GetElementsFromEdge(segGeom);
1056 "The periodic boundaries belong to "
1057 "more than one element of the mesh");
1060 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1063 allEdges[c->m_geomVec[i]->GetGlobalID()] =
1064 make_pair(elmt->at(0).second,
1065 geom->GetEorient(elmt->at(0).second));
1069 if (vComm->GetSize() == 1)
1071 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1074 vector<int> vertList(2);
1075 vertList[0] = segGeom->GetVid(0);
1076 vertList[1] = segGeom->GetVid(1);
1077 allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1080 if (vComm->GetSize() == 1)
1082 compOrder[it.second->begin()->first] = tmpOrder;
1087 if (perComps.count(cId1) == 0)
1089 if (perComps.count(cId2) == 0)
1091 perComps[cId1] = cId2;
1095 std::stringstream ss;
1096 ss <<
"Boundary region " << cId2 <<
" should be "
1097 <<
"periodic with " << perComps[cId2] <<
" but "
1098 <<
"found " << cId1 <<
" instead!";
1099 ASSERTL0(perComps[cId2] == cId1, ss.str());
1104 std::stringstream ss;
1105 ss <<
"Boundary region " << cId1 <<
" should be "
1106 <<
"periodic with " << perComps[cId1] <<
" but "
1107 <<
"found " << cId2 <<
" instead!";
1108 ASSERTL0(perComps[cId1] == cId1, ss.str());
1115 for (
auto &cIt : perComps)
1117 idmax = max(idmax, cIt.first);
1118 idmax = max(idmax, cIt.second);
1123 for (
auto &cIt : perComps)
1125 perid[cIt.first] = cIt.second;
1129 for (
int i = 0; i < idmax; ++i)
1135 if (perComps.count(perid[i]))
1139 perComps[i] = perid[i];
1145 int n = vComm->GetSize();
1146 int p = vComm->GetRank();
1152 edgecounts[
p] = allEdges.size();
1156 for (i = 1; i < n; ++i)
1158 edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
1167 auto sIt = allEdges.begin();
1169 for (i = 0; sIt != allEdges.end(); ++sIt)
1171 edgeIds[edgeoffset[
p] + i] = sIt->first;
1172 edgeIdx[edgeoffset[
p] + i] = sIt->second.first;
1173 edgeOrient[edgeoffset[
p] + i] = sIt->second.second;
1174 edgeVerts[edgeoffset[
p] + i++] = allVerts[sIt->first].size();
1197 for (i = 0; i < n; ++i)
1199 if (edgecounts[i] > 0)
1202 edgeVerts + edgeoffset[i], 1);
1211 for (i = 1; i < n; ++i)
1213 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1217 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1219 for (j = 0; j < allVerts[sIt->first].size(); ++j)
1221 traceIds[vertoffset[
p] + i++] = allVerts[sIt->first][j];
1228 map<int, StdRegions::Orientation> orientMap;
1229 map<int, vector<int>> vertMap;
1231 for (cnt = i = 0; i < totEdges; ++i)
1233 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1234 "Already found edge in orientation map!");
1254 orientMap[edgeIds[i]] = o;
1256 vector<int> verts(edgeVerts[i]);
1258 for (j = 0; j < edgeVerts[i]; ++j)
1260 verts[j] = traceIds[cnt++];
1262 vertMap[edgeIds[i]] = verts;
1269 map<int, int> allCompPairs;
1278 for (
auto &cIt : perComps)
1281 const int id1 = cIt.first;
1282 const int id2 = cIt.second;
1283 std::string id1s = std::to_string(id1);
1284 std::string id2s = std::to_string(id2);
1286 if (compMap.count(id1) > 0)
1288 c[0] = compMap[id1];
1291 if (compMap.count(id2) > 0)
1293 c[1] = compMap[id2];
1299 map<int, int> compPairs;
1302 "Unable to find composite " + id1s +
" in order map.");
1304 "Unable to find composite " + id2s +
" in order map.");
1305 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1306 "Periodic composites " + id1s +
" and " + id2s +
1307 " should have the same number of elements.");
1308 ASSERTL0(compOrder[id1].size() > 0,
"Periodic composites " +
1309 id1s +
" and " + id2s +
1313 for (i = 0; i < compOrder[id1].size(); ++i)
1315 int eId1 = compOrder[id1][i];
1316 int eId2 = compOrder[id2][i];
1318 ASSERTL0(compPairs.count(eId1) == 0,
"Already paired.");
1320 if (compPairs.count(eId2) != 0)
1322 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1324 compPairs[eId1] = eId2;
1330 for (i = 0; i < 2; ++i)
1337 if (c[i]->m_geomVec.size() > 0)
1339 for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1341 locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1347 for (
auto &pIt : compPairs)
1349 int ids[2] = {pIt.first, pIt.second};
1350 bool local[2] = {locEdges.count(pIt.first) > 0,
1351 locEdges.count(pIt.second) > 0};
1353 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1354 orientMap.count(ids[1]) > 0,
1355 "Unable to find edge in orientation map.");
1357 allCompPairs[pIt.first] = pIt.second;
1358 allCompPairs[pIt.second] = pIt.first;
1360 for (i = 0; i < 2; ++i)
1367 int other = (i + 1) % 2;
1370 orientMap[ids[i]] == orientMap[ids[other]]
1378 for (i = 0; i < 2; ++i)
1380 int other = (i + 1) % 2;
1383 orientMap[ids[i]] == orientMap[ids[other]]
1388 vector<int> perVerts1 = vertMap[ids[i]];
1389 vector<int> perVerts2 = vertMap[ids[other]];
1391 map<int, pair<int, bool>> tmpMap;
1394 tmpMap[perVerts1[0]] = make_pair(
1395 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1396 tmpMap[perVerts1[1]] = make_pair(
1397 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1401 tmpMap[perVerts1[0]] = make_pair(
1402 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1403 tmpMap[perVerts1[1]] = make_pair(
1404 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1407 for (
auto &mIt : tmpMap)
1413 auto perIt = periodicVerts.find(mIt.first);
1415 if (perIt == periodicVerts.end())
1417 periodicVerts[mIt.first].push_back(ent2);
1418 perIt = periodicVerts.find(mIt.first);
1423 for (j = 0; j < perIt->second.size(); ++j)
1425 if (perIt->second[j].id == mIt.second.first)
1434 perIt->second.push_back(ent2);
1445 for (cnt = i = 0; i < totEdges; ++i)
1447 int edgeId = edgeIds[i];
1449 ASSERTL0(allCompPairs.count(edgeId) > 0,
1450 "Unable to find matching periodic edge.");
1452 int perEdgeId = allCompPairs[edgeId];
1454 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1456 int vId = traceIds[cnt];
1458 auto perId = periodicVerts.find(vId);
1460 if (perId == periodicVerts.end())
1469 orientMap[edgeId] == orientMap[perEdgeId]
1470 ? vertMap[perEdgeId][(j + 1) % 2]
1471 : vertMap[perEdgeId][j];
1475 locVerts.count(perVertexId) > 0);
1477 periodicVerts[vId].push_back(ent);
1484 for (
auto &perIt : periodicVerts)
1487 for (i = 0; i < perIt.second.size(); ++i)
1489 auto perIt2 = periodicVerts.find(perIt.second[i].id);
1490 ASSERTL0(perIt2 != periodicVerts.end(),
1491 "Couldn't find periodic vertex.");
1493 for (j = 0; j < perIt2->second.size(); ++j)
1495 if (perIt2->second[j].id == perIt.first)
1501 for (k = 0; k < perIt.second.size(); ++k)
1503 if (perIt2->second[j].id == perIt.second[k].id)
1512 perIt.second.push_back(perIt2->second[j]);
1520 for (
auto &perIt : periodicVerts)
1522 if (locVerts.count(perIt.first) > 0)
1532 m_graph->GetCompositeOrdering();
1534 m_graph->GetBndRegionOrdering();
1549 map<int, RotPeriodicInfo> rotComp;
1550 map<int, int> perComps;
1551 map<int, vector<int>> allVerts;
1552 map<int, SpatialDomains::PointGeomVector> allCoord;
1553 map<int, vector<int>> allEdges;
1554 map<int, vector<StdRegions::Orientation>> allOrient;
1559 int region1ID, region2ID, i, j, k, cnt;
1563 for (i = 0; i < (*m_exp).size(); ++i)
1565 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1567 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1568 locVerts.insert(
id);
1571 for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1573 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1574 locEdges.insert(
id);
1581 for (
auto &it : bregions)
1587 if (locBCond->GetBoundaryConditionType() !=
1594 region1ID = it.first;
1596 std::static_pointer_cast<
1598 ->m_connectedBoundaryRegion;
1602 "Boundary region " + std::to_string(region1ID) +
1603 " should only contain 1 composite.");
1611 if (vComm->GetSize() == 1)
1613 cId1 = it.second->begin()->first;
1614 cId2 = bregions.find(region2ID)->second->begin()->first;
1618 cId1 = bndRegOrder.find(region1ID)->second[0];
1619 cId2 = bndRegOrder.find(region2ID)->second[0];
1623 if (boost::icontains(locBCond->GetUserDefined(),
"Rotated"))
1625 vector<string> tmpstr;
1627 boost::split(tmpstr, locBCond->GetUserDefined(),
1628 boost::is_any_of(
":"));
1630 if (boost::iequals(tmpstr[0],
"Rotated"))
1633 "Expected Rotated user defined string to "
1634 "contain direction and rotation angle "
1635 "and optionally a tolerance, "
1636 "i.e. Rotated:dir:PI/2:1e-6");
1638 ASSERTL1((tmpstr[1] ==
"x") || (tmpstr[1] ==
"y") ||
1641 "not specified as x,y or z");
1644 RotInfo.
m_dir = (tmpstr[1] ==
"x") ? 0
1645 : (tmpstr[1] ==
"y") ? 1
1652 if (tmpstr.size() == 4)
1656 RotInfo.
m_tol = std::stod(tmpstr[3]);
1661 "failed to cast tolerance input "
1662 "to a double value in Rotated"
1663 "boundary information");
1668 RotInfo.
m_tol = 1e-8;
1670 rotComp[cId1] = RotInfo;
1675 it.second->begin()->second;
1677 vector<unsigned int> tmpOrder;
1685 for (i = 0; i < c->m_geomVec.size(); ++i)
1688 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1690 ASSERTL1(faceGeom,
"Unable to cast to shared ptr");
1693 int faceId = c->m_geomVec[i]->GetGlobalID();
1694 locFaces.insert(faceId);
1698 if (vComm->GetSize() == 1)
1700 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1705 vector<int> vertList, edgeList;
1707 vector<StdRegions::Orientation> orientVec;
1708 for (j = 0; j < faceGeom->GetNumVerts(); ++j)
1710 vertList.push_back(faceGeom->GetVid(j));
1711 edgeList.push_back(faceGeom->GetEid(j));
1712 coordVec.push_back(faceGeom->GetVertex(j));
1713 orientVec.push_back(faceGeom->GetEorient(j));
1716 allVerts[faceId] = vertList;
1717 allEdges[faceId] = edgeList;
1718 allCoord[faceId] = coordVec;
1719 allOrient[faceId] = orientVec;
1724 if (vComm->GetSize() == 1)
1726 compOrder[it.second->begin()->first] = tmpOrder;
1732 if (perComps.count(cId1) == 0)
1734 if (perComps.count(cId2) == 0)
1736 perComps[cId1] = cId2;
1740 std::stringstream ss;
1741 ss <<
"Boundary region " << cId2 <<
" should be "
1742 <<
"periodic with " << perComps[cId2] <<
" but "
1743 <<
"found " << cId1 <<
" instead!";
1744 ASSERTL0(perComps[cId2] == cId1, ss.str());
1749 std::stringstream ss;
1750 ss <<
"Boundary region " << cId1 <<
" should be "
1751 <<
"periodic with " << perComps[cId1] <<
" but "
1752 <<
"found " << cId2 <<
" instead!";
1753 ASSERTL0(perComps[cId1] == cId1, ss.str());
1760 for (
auto &cIt : perComps)
1762 idmax = max(idmax, cIt.first);
1763 idmax = max(idmax, cIt.second);
1768 for (
auto &cIt : perComps)
1770 perid[cIt.first] = cIt.second;
1774 for (
int i = 0; i < idmax; ++i)
1780 if (perComps.count(perid[i]))
1784 perComps[i] = perid[i];
1791 int n = vComm->GetSize();
1792 int p = vComm->GetRank();
1802 rotcounts[
p] = rotComp.size();
1808 for (i = 1; i < n; ++i)
1810 rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
1819 auto rIt = rotComp.begin();
1821 for (i = 0; rIt != rotComp.end(); ++rIt)
1823 compid[rotoffset[
p] + i] = rIt->first;
1824 rotdir[rotoffset[
p] + i] = rIt->second.m_dir;
1825 rotangle[rotoffset[
p] + i] = rIt->second.m_angle;
1826 rottol[rotoffset[
p] + i++] = rIt->second.m_tol;
1835 for (i = 0; i < totrot; ++i)
1839 rotComp[compid[i]] = rinfo;
1844 facecounts[
p] = locFaces.size();
1850 for (i = 1; i < n; ++i)
1852 faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
1866 auto sIt = locFaces.begin();
1867 for (i = 0; sIt != locFaces.end(); ++sIt)
1869 faceIds[faceoffset[
p] + i] = *sIt;
1870 faceVerts[faceoffset[
p] + i++] = allVerts[*sIt].size();
1893 for (i = 0; i < n; ++i)
1895 if (facecounts[i] > 0)
1898 faceVerts + faceoffset[i], 1);
1909 for (i = 1; i < n; ++i)
1911 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1925 for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
1927 for (j = 0; j < allVerts[*sIt].size(); ++j)
1929 int vertId = allVerts[*sIt][j];
1930 vertIds[vertoffset[
p] + cnt] = vertId;
1931 vertX[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(0);
1932 vertY[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(1);
1933 vertZ[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(2);
1934 edgeIds[vertoffset[
p] + cnt] = allEdges[*sIt][j];
1935 edgeOrt[vertoffset[
p] + cnt++] = allOrient[*sIt][j];
1950 map<int, vector<int>> vertMap;
1951 map<int, vector<int>> edgeMap;
1952 map<int, SpatialDomains::PointGeomVector> coordMap;
1958 map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
1959 map<int, pair<int, int>> eIdMap;
1961 for (cnt = i = 0; i < totFaces; ++i)
1963 vector<int> edges(faceVerts[i]);
1964 vector<int> verts(faceVerts[i]);
1970 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1972 edges[j] = edgeIds[cnt];
1973 verts[j] = vertIds[cnt];
1977 vCoMap[vertIds[cnt]] = coord[j];
1980 auto testIns = eIdMap.insert(make_pair(
1982 make_pair(vertIds[tmp + j],
1983 vertIds[tmp + ((j + 1) % faceVerts[i])])));
1985 if (testIns.second ==
false)
2008 swap(testIns.first->second.first,
2009 testIns.first->second.second);
2013 vertMap[faceIds[i]] = verts;
2014 edgeMap[faceIds[i]] = edges;
2015 coordMap[faceIds[i]] = coord;
2033 map<int, map<StdRegions::Orientation, vector<int>>> vmap;
2034 map<int, map<StdRegions::Orientation, vector<int>>> emap;
2036 map<StdRegions::Orientation, vector<int>> quadVertMap;
2046 map<StdRegions::Orientation, vector<int>> quadEdgeMap;
2056 map<StdRegions::Orientation, vector<int>> triVertMap;
2060 map<StdRegions::Orientation, vector<int>> triEdgeMap;
2064 vmap[3] = triVertMap;
2065 vmap[4] = quadVertMap;
2066 emap[3] = triEdgeMap;
2067 emap[4] = quadEdgeMap;
2069 map<int, int> allCompPairs;
2073 map<int, int> fIdToCompId;
2078 for (
auto &cIt : perComps)
2081 const int id1 = cIt.first;
2082 const int id2 = cIt.second;
2083 std::string id1s = std::to_string(id1);
2084 std::string id2s = std::to_string(id2);
2086 if (compMap.count(id1) > 0)
2088 c[0] = compMap[id1];
2091 if (compMap.count(id2) > 0)
2093 c[1] = compMap[id2];
2099 map<int, int> compPairs;
2102 "Unable to find composite " + id1s +
" in order map.");
2104 "Unable to find composite " + id2s +
" in order map.");
2105 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
2106 "Periodic composites " + id1s +
" and " + id2s +
2107 " should have the same number of elements.");
2108 ASSERTL0(compOrder[id1].size() > 0,
"Periodic composites " +
2109 id1s +
" and " + id2s +
2113 for (i = 0; i < compOrder[id1].size(); ++i)
2115 int eId1 = compOrder[id1][i];
2116 int eId2 = compOrder[id2][i];
2118 ASSERTL0(compPairs.count(eId1) == 0,
"Already paired.");
2121 if (compPairs.count(eId2) != 0)
2123 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
2125 compPairs[eId1] = eId2;
2128 fIdToCompId[eId1] = id1;
2129 fIdToCompId[eId2] = id2;
2135 for (
auto &pIt : compPairs)
2137 int ids[2] = {pIt.first, pIt.second};
2138 bool local[2] = {locFaces.count(pIt.first) > 0,
2139 locFaces.count(pIt.second) > 0};
2141 ASSERTL0(coordMap.count(ids[0]) > 0 &&
2142 coordMap.count(ids[1]) > 0,
2143 "Unable to find face in coordinate map");
2145 allCompPairs[pIt.first] = pIt.second;
2146 allCompPairs[pIt.second] = pIt.first;
2151 coordMap[ids[0]], coordMap[ids[1]]};
2153 ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
2154 "Two periodic faces have different number "
2163 bool rotbnd =
false;
2170 if (rotComp.count(fIdToCompId[pIt.first]))
2173 dir = rotComp[fIdToCompId[pIt.first]].m_dir;
2174 angle = rotComp[fIdToCompId[pIt.first]].m_angle;
2175 tol = rotComp[fIdToCompId[pIt.first]].m_tol;
2179 for (i = 0; i < 2; ++i)
2187 int other = (i + 1) % 2;
2190 sign = (i == 0) ? 1.0 : -1.0;
2193 if (tmpVec[0].size() == 3)
2196 tmpVec[i], tmpVec[other], rotbnd, dir,
2202 tmpVec[i], tmpVec[other], rotbnd, dir,
2212 int nFaceVerts = vertMap[ids[0]].size();
2215 for (i = 0; i < 2; ++i)
2217 int other = (i + 1) % 2;
2220 sign = (i == 0) ? 1.0 : -1.0;
2223 if (tmpVec[0].size() == 3)
2226 tmpVec[i], tmpVec[other], rotbnd, dir,
2232 tmpVec[i], tmpVec[other], rotbnd, dir,
2236 if (nFaceVerts == 3)
2241 "Unsupported face orientation for face " +
2242 std::to_string(ids[i]));
2246 vector<int> per1 = vertMap[ids[i]];
2247 vector<int> per2 = vertMap[ids[other]];
2251 map<int, pair<int, bool>> tmpMap;
2255 for (j = 0; j < nFaceVerts; ++j)
2257 int v = vmap[nFaceVerts][o][j];
2259 make_pair(per2[v], locVerts.count(per2[v]) > 0);
2263 for (
auto &mIt : tmpMap)
2270 auto perIt = periodicVerts.find(mIt.first);
2272 if (perIt == periodicVerts.end())
2276 periodicVerts[mIt.first].push_back(ent2);
2277 perIt = periodicVerts.find(mIt.first);
2284 for (k = 0; k < perIt->second.size(); ++k)
2286 if (perIt->second[k].id == mIt.second.first)
2292 if (k == perIt->second.size())
2294 perIt->second.push_back(ent2);
2302 for (i = 0; i < 2; ++i)
2304 int other = (i + 1) % 2;
2307 sign = (i == 0) ? 1.0 : -1.0;
2309 if (tmpVec[0].size() == 3)
2312 tmpVec[i], tmpVec[other], rotbnd, dir,
2318 tmpVec[i], tmpVec[other], rotbnd, dir,
2322 vector<int> per1 = edgeMap[ids[i]];
2323 vector<int> per2 = edgeMap[ids[other]];
2325 map<int, pair<int, bool>> tmpMap;
2327 for (j = 0; j < nFaceVerts; ++j)
2329 int e = emap[nFaceVerts][o][j];
2331 make_pair(per2[e], locEdges.count(per2[e]) > 0);
2334 for (
auto &mIt : tmpMap)
2341 auto perIt = periodicEdges.find(mIt.first);
2343 if (perIt == periodicEdges.end())
2345 periodicEdges[mIt.first].push_back(ent2);
2346 perIt = periodicEdges.find(mIt.first);
2350 for (k = 0; k < perIt->second.size(); ++k)
2352 if (perIt->second[k].id == mIt.second.first)
2358 if (k == perIt->second.size())
2360 perIt->second.push_back(ent2);
2370 ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
2371 "At this point the size of allCompPairs "
2372 "should have been the same as fIdToCompId");
2376 map<int, int> eIdToCompId;
2382 for (cnt = i = 0; i < totFaces; ++i)
2384 bool rotbnd =
false;
2389 int faceId = faceIds[i];
2391 ASSERTL0(allCompPairs.count(faceId) > 0,
2392 "Unable to find matching periodic face. faceId = " +
2393 std::to_string(faceId));
2395 int perFaceId = allCompPairs[faceId];
2398 ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
2399 "Face " + std::to_string(faceId) +
2400 " not found in fIdtoCompId map");
2401 if (rotComp.count(fIdToCompId[faceId]))
2404 dir = rotComp[fIdToCompId[faceId]].m_dir;
2405 angle = rotComp[fIdToCompId[faceId]].m_angle;
2406 tol = rotComp[fIdToCompId[faceId]].m_tol;
2409 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
2411 int vId = vertIds[cnt];
2413 auto perId = periodicVerts.find(vId);
2415 if (perId == periodicVerts.end())
2424 coordMap[faceId], coordMap[perFaceId]};
2426 int nFaceVerts = tmpVec[0].size();
2430 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2433 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2439 vertMap[perFaceId][vmap[nFaceVerts][o][j]];
2443 locVerts.count(perVertexId) > 0);
2445 periodicVerts[vId].push_back(ent);
2448 int eId = edgeIds[cnt];
2450 perId = periodicEdges.find(eId);
2456 eIdToCompId[eId] = fIdToCompId[faceId];
2459 if (perId == periodicEdges.end())
2467 coordMap[faceId], coordMap[perFaceId]};
2469 int nFaceEdges = tmpVec[0].size();
2473 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2476 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2482 edgeMap[perFaceId][emap[nFaceEdges][o][j]];
2485 locEdges.count(perEdgeId) > 0);
2487 periodicEdges[eId].push_back(ent);
2493 eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
2501 for (
auto &perIt : periodicVerts)
2504 for (i = 0; i < perIt.second.size(); ++i)
2507 auto perIt2 = periodicVerts.find(perIt.second[i].id);
2508 ASSERTL0(perIt2 != periodicVerts.end(),
2509 "Couldn't find periodic vertex.");
2514 for (j = 0; j < perIt2->second.size(); ++j)
2516 if (perIt2->second[j].id == perIt.first)
2521 for (k = 0; k < perIt.second.size(); ++k)
2523 if (perIt2->second[j].id == perIt.second[k].id)
2529 if (k == perIt.second.size())
2531 perIt.second.push_back(perIt2->second[j]);
2537 for (
auto &perIt : periodicEdges)
2539 for (i = 0; i < perIt.second.size(); ++i)
2541 auto perIt2 = periodicEdges.find(perIt.second[i].id);
2542 ASSERTL0(perIt2 != periodicEdges.end(),
2543 "Couldn't find periodic edge.");
2545 for (j = 0; j < perIt2->second.size(); ++j)
2547 if (perIt2->second[j].id == perIt.first)
2552 for (k = 0; k < perIt.second.size(); ++k)
2554 if (perIt2->second[j].id == perIt.second[k].id)
2560 if (k == perIt.second.size())
2562 perIt.second.push_back(perIt2->second[j]);
2569 for (
auto &perIt : periodicEdges)
2571 bool rotbnd =
false;
2577 auto eIt = eIdMap.find(perIt.first);
2579 *vCoMap[eIt->second.second]};
2582 if (rotComp.count(eIdToCompId[perIt.first]))
2585 dir = rotComp[eIdToCompId[perIt.first]].m_dir;
2586 angle = rotComp[eIdToCompId[perIt.first]].m_angle;
2587 tol = rotComp[eIdToCompId[perIt.first]].m_tol;
2593 for (i = 0; i < perIt.second.size(); ++i)
2595 eIt = eIdMap.find(perIt.second[i].id);
2598 *vCoMap[eIt->second.first],
2599 *vCoMap[eIt->second.second]};
2601 int vMap[2] = {-1, -1};
2607 r.
Rotate(v[0], dir, angle);
2609 if (r.
dist(
w[0]) < tol)
2615 r.
Rotate(v[1], dir, angle);
2616 if (r.
dist(
w[0]) < tol)
2623 "Unable to align rotationally "
2624 "periodic edge vertex");
2631 0.5 * (
w[0](0) - v[0](0) +
w[1](0) - v[1](0));
2633 0.5 * (
w[0](1) - v[0](1) +
w[1](1) - v[1](1));
2635 0.5 * (
w[0](2) - v[0](2) +
w[1](2) - v[1](2));
2637 for (j = 0; j < 2; ++j)
2642 for (k = 0; k < 2; ++k)
2648 if (
sqrt((x1 - x) * (x1 - x) +
2649 (y1 - y) * (y1 - y) +
2650 (z1 -
z) * (z1 -
z)) < 1e-8)
2659 ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
2660 "Unable to align periodic edge vertex.");
2661 ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
2662 (vMap[1] == 0 || vMap[1] == 1) &&
2663 (vMap[0] != vMap[1]),
2664 "Unable to align periodic edge vertex.");
2678 for (
auto &perIt : periodicVerts)
2680 if (locVerts.count(perIt.first) > 0)
2686 for (
auto &perIt : periodicEdges)
2688 if (locEdges.count(perIt.first) > 0)
2696 ASSERTL1(
false,
"Not setup for this expansion");
2708 "Routine currently only tested for HybridDGHelmholtz");
2711 "Full matrix global systems are not supported for HDG "
2715 "The local to global map is not set up for the requested "
2724 (*m_globalBndMat)[mkey] = glo_matrix;
2728 glo_matrix = matrixIter->second;
2743 bool returnval =
true;
2762 int vSame = (returnval ? 1 : 0);
2765 return (vSame == 1);
2780 for (
int v = 0; v < 2; ++v)
2785 if (vertExp->GetLeftAdjacentElementExp()
2788 (*m_exp)[i]->GetGeom()->GetGlobalID())
2920 "Routine not set up for Gauss Lagrange points distribution");
2937 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2949 bool PutFwdInBwdOnBCs,
bool DoExchange)
2981 for (
int n = 0, cnt = 0; n < nexp; ++n)
2986 for (
int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2989 m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2994 exp->GetTracePhysVals(e, elmtToTrace[n][e],
field + phys_offset,
3010 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
3021 bool PutFwdInBwdOnBCs)
3034 bool PutFwdInBwdOnBCs)
3038 if (PutFwdInBwdOnBCs)
3041 for (
int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3043 if (bndConditions[n]->GetBoundaryConditionType() ==
3046 auto ne = bndCondExpansions[n]->GetExpSize();
3047 for (
int e = 0; e < ne; ++e)
3049 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3050 auto id2 =
m_trace->GetPhys_Offset(
3051 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3057 else if (bndConditions[n]->GetBoundaryConditionType() ==
3059 bndConditions[n]->GetBoundaryConditionType() ==
3062 auto ne = bndCondExpansions[n]->GetExpSize();
3063 for (
int e = 0; e < ne; ++e)
3065 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3066 auto id2 =
m_trace->GetPhys_Offset(
3067 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3072 else if (bndConditions[n]->GetBoundaryConditionType() !=
3075 ASSERTL0(
false,
"Method not set up for this "
3076 "boundary condition.");
3082 for (
int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3084 if (bndConditions[n]->GetBoundaryConditionType() ==
3087 auto ne = bndCondExpansions[n]->GetExpSize();
3088 for (
int e = 0; e < ne; ++e)
3090 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3091 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3092 auto id2 =
m_trace->GetPhys_Offset(
3093 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3100 else if (bndConditions[n]->GetBoundaryConditionType() ==
3102 bndConditions[n]->GetBoundaryConditionType() ==
3106 for (
int e = 0; e < ne; ++e)
3108 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3109 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3112 "Method not set up for non-zero "
3113 "Neumann boundary condition");
3114 auto id2 =
m_trace->GetPhys_Offset(
3115 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3122 else if (bndConditions[n]->GetBoundaryConditionType() ==
3127 else if (bndConditions[n]->GetBoundaryConditionType() !=
3131 "Method not set up for this boundary "
3171 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3206 m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3212 int n,
p, offset, phys_offset;
3217 "input array is of insufficient length");
3218 for (n = 0; n < nexp; ++n)
3221 for (
p = 0;
p < (*m_exp)[n]->GetNtraces(); ++
p)
3224 m_trace->GetPhys_Offset(elmtToTrace[n][
p]->GetElmtId());
3225 (*m_exp)[n]->GetTracePhysVals(
p, elmtToTrace[n][
p],
3226 inarray + phys_offset,
3227 t_tmp = outarray + offset);
3255 int n, offset, t_offset;
3262 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3265 if ((*
m_exp)[n]->GetBasis(0)->GetBasisType() !=
3269 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3270 if (negatedFluxNormal[2 * n])
3272 outarray[offset] -= Fn[t_offset];
3276 outarray[offset] += Fn[t_offset];
3279 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3280 if (negatedFluxNormal[2 * n + 1])
3282 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3287 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3295 static int sav_ncoeffs = 0;
3296 if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3306 m_Ixm = BASE->GetI(coords);
3308 m_Ixp = BASE->GetI(coords);
3309 sav_ncoeffs = e_ncoeffs;
3312 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3313 if (negatedFluxNormal[2 * n])
3315 for (j = 0; j < e_ncoeffs; j++)
3317 outarray[offset + j] -=
3318 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3323 for (j = 0; j < e_ncoeffs; j++)
3325 outarray[offset + j] +=
3326 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3330 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3331 if (negatedFluxNormal[2 * n + 1])
3333 for (j = 0; j < e_ncoeffs; j++)
3335 outarray[offset + j] -=
3336 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3341 for (j = 0; j < e_ncoeffs; j++)
3343 outarray[offset + j] +=
3344 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3358 m_trace->IProductWRTBase(Fn, Fcoeffs);
3364 int e, n, offset, t_offset;
3373 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3375 t_offset =
GetTrace()->GetPhys_Offset(
3376 elmtToTrace[n][e]->GetElmtId());
3377 (*m_exp)[n]->AddEdgeNormBoundaryInt(
3378 e, elmtToTrace[n][e], Fn + t_offset,
3379 e_outarray = outarray + offset);
3388 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3390 t_offset =
GetTrace()->GetPhys_Offset(
3391 elmtToTrace[n][e]->GetElmtId());
3392 (*m_exp)[n]->AddFaceNormBoundaryInt(
3393 e, elmtToTrace[n][e], Fn + t_offset,
3394 e_outarray = outarray + offset);
3426 "tested for 1D expansion");
3428 m_trace->IProductWRTBase(Fwd, Coeffs);
3430 m_trace->IProductWRTBase(Bwd, Coeffs);
3440 const bool PhysSpaceForcing)
3442 int i, n, cnt, nbndry;
3450 if (PhysSpaceForcing)
3462 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
3463 int NumDirBCs =
m_traceMap->GetNumLocalDirBndCoeffs();
3470 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
3479 for (cnt = n = 0; n < nexp; ++n)
3481 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3482 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3488 Floc =
Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3492 m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3508 locid = bndCondMap[cnt + j];
3509 loclambda[locid] = Sign[locid] * bndcoeffs[j];
3520 locid = bndCondMap[cnt + j];
3521 bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3538 if (GloBndDofs - NumDirBCs > 0)
3543 LinSys->Solve(bndrhs, loclambda,
m_traceMap);
3561 out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3603 const std::string varName,
3619 npoints = locExpList->GetNpoints();
3624 locExpList->GetCoords(x0, x1, x2);
3648 0, (std::static_pointer_cast<
3651 ->m_dirichletCondition)
3652 .Evaluate(x0[0], x1[0], x2[0], time));
3660 0, (std::static_pointer_cast<
3663 ->m_neumannCondition)
3664 .Evaluate(x0[0], x1[0], x2[0], time));
3670 0, (std::static_pointer_cast<
3674 .Evaluate(x0[0], x1[0], x2[0], time));
3688 "This type of BC not implemented yet");
3697 std::static_pointer_cast<
3702 valuesExp(npoints, 1.0);
3704 string bcfilename = bcPtr->m_filename;
3705 string exprbcs = bcPtr->m_expr;
3707 if (bcfilename !=
"")
3709 locExpList->ExtractCoeffsFromFile(
3710 bcfilename, bcPtr->GetComm(), varName,
3711 locExpList->UpdateCoeffs());
3712 locExpList->BwdTrans(locExpList->GetCoeffs(),
3713 locExpList->UpdatePhys());
3714 valuesFile = locExpList->GetPhys();
3720 std::static_pointer_cast<
3723 ->m_dirichletCondition;
3725 condition.
Evaluate(x0, x1, x2, time, valuesExp);
3729 locExpList->UpdatePhys(), 1);
3730 locExpList->FwdTransBndConstrained(
3731 locExpList->GetPhys(), locExpList->UpdateCoeffs());
3737 std::static_pointer_cast<
3742 if (bcfilename !=
"")
3744 locExpList->ExtractCoeffsFromFile(
3745 bcfilename, bcPtr->GetComm(), varName,
3746 locExpList->UpdateCoeffs());
3747 locExpList->BwdTrans(locExpList->GetCoeffs(),
3748 locExpList->UpdatePhys());
3753 std::static_pointer_cast<
3756 ->m_neumannCondition;
3757 condition.
Evaluate(x0, x1, x2, time,
3758 locExpList->UpdatePhys());
3761 locExpList->IProductWRTBase(locExpList->GetPhys(),
3762 locExpList->UpdateCoeffs());
3768 std::static_pointer_cast<
3774 if (bcfilename !=
"")
3776 locExpList->ExtractCoeffsFromFile(
3777 bcfilename, bcPtr->GetComm(), varName,
3778 locExpList->UpdateCoeffs());
3779 locExpList->BwdTrans(locExpList->GetCoeffs(),
3780 locExpList->UpdatePhys());
3785 std::static_pointer_cast<
3789 condition.
Evaluate(x0, x1, x2, time,
3790 locExpList->UpdatePhys());
3793 locExpList->IProductWRTBase(locExpList->GetPhys(),
3794 locExpList->UpdateCoeffs());
3804 "This type of BC not implemented yet");
3833 id2 =
m_trace->GetPhys_Offset(
3834 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3849 id2 =
m_trace->GetPhys_Offset(
3850 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3861 "Method not set up for this boundary condition.");
3878 map<int, int> VertGID;
3903 VertGID[Vid] = cnt++;
3912 for (i = 0; i < exp->GetNverts(); ++i)
3914 id = exp->GetGeom()->GetVid(i);
3916 if (VertGID.count(
id) > 0)
3918 bid = VertGID.find(
id)->second;
3920 "Edge already set");
3928 "Failed to visit all boundary condtiions");
3933 map<int, int> globalIdMap;
3942 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3968 m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3971 globalIdMap[(*tmp)[0].first->GetGlobalID()];
3979 map<int, int> globalIdMap;
3989 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
4013 m_graph->GetElementsFromFace(exp2d->GetGeom2D());
4015 globalIdMap[tmp->at(0).first->GetGlobalID()];
4022 ASSERTL1(
false,
"Not setup for this expansion");
4032 std::shared_ptr<ExpList> &result,
4033 const bool DeclareCoeffPhysArrays)
4036 std::vector<unsigned int> eIDs;
4042 for (cnt = n = 0; n < i; ++n)
4050 eIDs.push_back(ElmtID[cnt + n]);
4057 if (DeclareCoeffPhysArrays)
4060 int offsetOld, offsetNew;
4062 for (n = 0; n < result->GetExpSize(); ++n)
4064 nq =
GetExp(ElmtID[cnt + n])->GetTotPoints();
4066 offsetNew = result->GetPhys_Offset(n);
4068 tmp2 = result->UpdatePhys() + offsetNew, 1);
4069 nq =
GetExp(ElmtID[cnt + n])->GetNcoeffs();
4071 offsetNew = result->GetCoeff_Offset(n);
4073 tmp2 = result->UpdateCoeffs() + offsetNew, 1);
4109 map<int, RobinBCInfoSharedPtr> returnval;
4125 int npoints = locExpList->GetNpoints();
4131 locExpList->GetCoords(x0, x1, x2);
4134 std::static_pointer_cast<
4136 ->m_robinPrimitiveCoeff;
4139 coeffeqn.
Evaluate(x0, x1, x2, 0.0, coeffphys);
4141 for (e = 0; e < locExpList->GetExpSize(); ++e)
4146 Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4148 elmtid = ElmtID[cnt + e];
4150 if (returnval.count(elmtid) != 0)
4152 rInfo->next = returnval.find(elmtid)->second;
4154 returnval[elmtid] = rInfo;
4177 int i, e, ncoeff_edge;
4187 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4191 edge_lambda = loc_lambda;
4199 int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4202 for (e = 0; e < nEdges; ++e)
4204 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4205 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4207 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4208 elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4210 edge_lambda = edge_lambda + ncoeff_edge;
4213 (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs +
m_coeff_offset[i],
4214 elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4215 cnt += (*m_exp)[i]->GetNcoeffs();
4248 int i, cnt, e, ncoeff_trace;
4255 int nq_elmt, nm_elmt;
4256 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4261 trace_lambda = loc_lambda;
4265 int num_points[] = {0, 0, 0};
4266 int num_modes[] = {0, 0, 0};
4271 nq_elmt = (*m_exp)[i]->GetTotPoints();
4272 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4276 out_tmp = force + nm_elmt;
4279 for (
int j = 0; j < dim; ++j)
4281 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4282 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4288 int nTraces = (*m_exp)[i]->GetNtraces();
4291 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4293 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4294 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4296 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4299 elmtToTrace[i][e]->SetCoeffsToOrientation(
4300 edgedir, traceCoeffs[e], traceCoeffs[e]);
4306 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4308 trace_lambda = trace_lambda + ncoeff_trace;
4323 num_modes[0], PkeyQ1);
4325 num_modes[1], PkeyQ2);
4327 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4328 (*
m_exp)[i]->GetGeom());
4330 BkeyQ1, BkeyQ2, qGeom);
4338 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4340 num_modes[0], PkeyT1);
4342 num_modes[1], PkeyT2);
4344 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4345 (*
m_exp)[i]->GetGeom());
4347 BkeyT1, BkeyT2, tGeom);
4359 num_modes[0], PkeyH1);
4361 num_modes[1], PkeyH2);
4363 num_modes[2], PkeyH3);
4365 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4366 (*
m_exp)[i]->GetGeom());
4368 BkeyH1, BkeyH2, BkeyH3, hGeom);
4376 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4378 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4380 num_modes[0], PkeyT1);
4382 num_modes[1], PkeyT2);
4384 num_modes[2], PkeyT3);
4386 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4387 (*
m_exp)[i]->GetGeom());
4389 BkeyT1, BkeyT2, BkeyT3, tGeom);
4394 "Wrong shape type, HDG postprocessing is not "
4401 elmtToTrace[i], traceCoeffs, out_tmp);
4402 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4403 ppExp->IProductWRTDerivBase(0, qrhs, force);
4407 elmtToTrace[i], traceCoeffs, out_tmp);
4409 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4410 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4412 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4415 (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs +
m_coeff_offset[i], qrhs);
4416 force[0] = (*m_exp)[i]->Integral(qrhs);
4421 ppExp->DetShapeType(), *ppExp);
4427 out = (*lapsys) * in;
4431 ppExp->BwdTrans(out.
GetPtr(), work);
4432 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray +
m_coeff_offset[i]);
4461 m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4463 m_trace->IProductWRTBase(BwdFlux, FCoeffs);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
#define sign(a, b)
return the sign(b)*a
Describes the specification for a Basis.
NekDouble Evaluate() const
Interpreter class for the evaluation of mathematical expressions.
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
Defines a specification for a set of points.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansi...
~DisContField() override
Destructor.
Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions() override
std::vector< int > m_periodicBwdCopy
void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
Add trace contributions into elemental coefficient spaces.
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
void SetUpDG(const std::string="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
Set up all DG member variables and maps.
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
std::set< int > m_boundaryTraces
A set storing the global IDs of any boundary Verts.
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
std::vector< bool > m_negatedFluxNormal
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Array< OneD, int > m_BCtoTraceMap
void EvaluateHDGPostProcessing(const Array< OneD, const NekDouble > &coeffs, Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces) override
Obtain a copy of the periodic edges and vertices for this field.
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions() override
MultiRegions::ExpListSharedPtr & v_UpdateBndCondExpansion(int i) override
void v_SetBndCondBwdWeight(const int index, const NekDouble value) override
void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID) override
void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray.
InterfaceMapDGSharedPtr & v_GetInterfaceMap(void) override
InterfaceMapDGSharedPtr m_interfaceMap
Interfaces mapping for trace space.
void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays) override
void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp) override
Fill the weight with m_bndCondBndWeight.
const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight() override
const Array< OneD, const MultiRegions::ExpListSharedPtr > & v_GetBndCondExpansions() override
void FindPeriodicTraces(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
Array< OneD, NekDouble > m_bndCondBndWeight
bool SameTypeOfBoundaryConditions(const DisContField &In)
Check to see if expansion has the same BCs as In.
ExpListSharedPtr & v_GetTrace() override
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition structure definition on the diff...
void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble) override
Evaluate all boundary conditions at a given time..
void FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions, const Array< OneD, const ExpListSharedPtr > &BndCondExpansions, bool PutFwdInBwdOnBCs)
void GetFwdBwdTracePhys(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndCond, const Array< OneD, const ExpListSharedPtr > &BndCondExp)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray) override
Add trace contributions into elemental coefficient spaces.
const LocTraceToTraceMapSharedPtr & v_GetLocTraceToTraceMap(void) const override
void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray) override
void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs) override
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
ExpListSharedPtr m_trace
Trace space storage for points between elements.
void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field) override
void v_Reset() override
Reset this field, so that geometry information can be updated.
std::vector< bool > & GetNegatedFluxNormal(void)
DisContField()
Default constructor.
std::vector< bool > & v_GetLeftAdjacentTraces(void) override
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
GlobalLinSysKey v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing) override
Solve the Helmholtz equation.
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
bool IsLeftAdjacentTrace(const int n, const int e)
This routine determines if an element is to the "left" of the adjacent trace, which arises from the i...
std::vector< bool > m_leftAdjacentTraces
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &coeffs, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd) override
Array< OneD, int > m_BCtoElmMap
AssemblyMapDGSharedPtr & v_GetTraceMap(void) override
std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo() override
void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true)
Discretises the boundary conditions.
Base class for all multi-elemental spectral/hp expansions.
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
std::shared_ptr< ExpList > & GetTrace()
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
int GetCoeff_Offset(int n) const
Get the start offset position for a local contiguous list of coeffs correspoinding to element n.
int GetExpSize(void)
This function returns the number of elements in the expansion.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
std::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...
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
bool m_physState
The state of the array m_phys.
LibUtilities::CommSharedPtr m_comm
Communicator.
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global This function calculates the error with respect to...
void ApplyGeomInfo()
Apply geometry information to each expansion.
LibUtilities::SessionReaderSharedPtr m_session
Session.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
int GetPhys_Offset(int n) const
Get the start offset position for a local contiguous list of quadrature points in a full array corres...
ExpansionType m_expType
Expansion type.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Describe a linear system.
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Array< OneD, DataType > & GetPtr()
const BoundaryRegionCollection & GetBoundaryRegions(void) const
const BoundaryConditionCollection & GetBoundaryConditions(void) const
void Rotate(PointGeom &a, int dir, NekDouble angle)
_this = rotation of a by angle 'angle' around axis dir
NekDouble dist(PointGeom &a)
return distance between this and input a
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2, bool doRot=false, int dir=0, NekDouble angle=0.0, NekDouble tol=1e-8)
Get the orientation of face1.
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2, bool doRot, int dir, NekDouble angle, NekDouble tol)
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
@ eOrtho_A
Principle Orthogonal Functions .
@ eOrtho_C
Principle Orthogonal Functions .
@ eOrtho_B
Principle Orthogonal Functions .
std::shared_ptr< Expansion > ExpansionSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
std::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
static ExpListSharedPtr NullExpListSharedPtr
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
std::shared_ptr< InterfaceMapDG > InterfaceMapDGSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static GlobalLinSysKey NullGlobalLinSysKey(StdRegions::eNoMatrixType)
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
static const NekDouble kNekUnsetDouble
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
std::map< int, std::vector< unsigned int > > CompositeOrdering
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
std::shared_ptr< std::vector< std::pair< GeometrySharedPtr, int > > > GeometryLinkSharedPtr
std::map< int, std::vector< unsigned int > > BndRegionOrdering
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
std::shared_ptr< Composite > CompositeSharedPtr
std::vector< PointGeomSharedPtr > PointGeomVector
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
std::shared_ptr< HexGeom > HexGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
std::map< int, CompositeSharedPtr > CompositeMap
@ eInvLaplacianWithUnityMean
std::map< ConstFactorType, NekDouble > ConstFactorMap
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1FwdDir1_Dir2FwdDir2
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
std::vector< double > w(NPUPPER)
std::vector< double > z(NPUPPER)
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
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.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
scalarT< T > sqrt(scalarT< T > in)
int id
Geometry ID of entity.
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
bool isLocal
Flag specifying if this entity is local to this partition.
NekDouble m_tol
Tolerance to rotation is considered identical.
int m_dir
Axis of rotation. 0 = 'x', 1 = 'y', 2 = 'z'.
NekDouble m_angle
Angle of rotation in radians.