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 " +
1037 boost::lexical_cast<string>(region1ID) +
1038 " should only contain 1 composite.");
1040 vector<unsigned int> tmpOrder;
1045 it.second->begin()->second;
1047 for (i = 0; i < c->m_geomVec.size(); ++i)
1050 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1052 ASSERTL0(segGeom,
"Unable to cast to shared ptr");
1055 m_graph->GetElementsFromEdge(segGeom);
1057 "The periodic boundaries belong to "
1058 "more than one element of the mesh");
1061 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1064 allEdges[c->m_geomVec[i]->GetGlobalID()] =
1065 make_pair(elmt->at(0).second,
1066 geom->GetEorient(elmt->at(0).second));
1070 if (vComm->GetSize() == 1)
1072 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1075 vector<int> vertList(2);
1076 vertList[0] = segGeom->GetVid(0);
1077 vertList[1] = segGeom->GetVid(1);
1078 allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1081 if (vComm->GetSize() == 1)
1083 compOrder[it.second->begin()->first] = tmpOrder;
1088 if (perComps.count(cId1) == 0)
1090 if (perComps.count(cId2) == 0)
1092 perComps[cId1] = cId2;
1096 std::stringstream ss;
1097 ss <<
"Boundary region " << cId2 <<
" should be "
1098 <<
"periodic with " << perComps[cId2] <<
" but "
1099 <<
"found " << cId1 <<
" instead!";
1100 ASSERTL0(perComps[cId2] == cId1, ss.str());
1105 std::stringstream ss;
1106 ss <<
"Boundary region " << cId1 <<
" should be "
1107 <<
"periodic with " << perComps[cId1] <<
" but "
1108 <<
"found " << cId2 <<
" instead!";
1109 ASSERTL0(perComps[cId1] == cId1, ss.str());
1116 for (
auto &cIt : perComps)
1118 idmax = max(idmax, cIt.first);
1119 idmax = max(idmax, cIt.second);
1124 for (
auto &cIt : perComps)
1126 perid[cIt.first] = cIt.second;
1130 for (
int i = 0; i < idmax; ++i)
1136 if (perComps.count(perid[i]))
1140 perComps[i] = perid[i];
1146 int n = vComm->GetSize();
1147 int p = vComm->GetRank();
1153 edgecounts[
p] = allEdges.size();
1157 for (i = 1; i < n; ++i)
1159 edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
1168 auto sIt = allEdges.begin();
1170 for (i = 0; sIt != allEdges.end(); ++sIt)
1172 edgeIds[edgeoffset[
p] + i] = sIt->first;
1173 edgeIdx[edgeoffset[
p] + i] = sIt->second.first;
1174 edgeOrient[edgeoffset[
p] + i] = sIt->second.second;
1175 edgeVerts[edgeoffset[
p] + i++] = allVerts[sIt->first].size();
1198 for (i = 0; i < n; ++i)
1200 if (edgecounts[i] > 0)
1203 edgeVerts + edgeoffset[i], 1);
1212 for (i = 1; i < n; ++i)
1214 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1218 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1220 for (j = 0; j < allVerts[sIt->first].size(); ++j)
1222 traceIds[vertoffset[
p] + i++] = allVerts[sIt->first][j];
1229 map<int, StdRegions::Orientation> orientMap;
1230 map<int, vector<int>> vertMap;
1232 for (cnt = i = 0; i < totEdges; ++i)
1234 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1235 "Already found edge in orientation map!");
1255 orientMap[edgeIds[i]] = o;
1257 vector<int> verts(edgeVerts[i]);
1259 for (j = 0; j < edgeVerts[i]; ++j)
1261 verts[j] = traceIds[cnt++];
1263 vertMap[edgeIds[i]] = verts;
1270 map<int, int> allCompPairs;
1279 for (
auto &cIt : perComps)
1282 const int id1 = cIt.first;
1283 const int id2 = cIt.second;
1284 std::string id1s = boost::lexical_cast<string>(id1);
1285 std::string id2s = boost::lexical_cast<string>(id2);
1287 if (compMap.count(id1) > 0)
1289 c[0] = compMap[id1];
1292 if (compMap.count(id2) > 0)
1294 c[1] = compMap[id2];
1300 map<int, int> compPairs;
1303 "Unable to find composite " + id1s +
" in order map.");
1305 "Unable to find composite " + id2s +
" in order map.");
1306 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1307 "Periodic composites " + id1s +
" and " + id2s +
1308 " should have the same number of elements.");
1309 ASSERTL0(compOrder[id1].size() > 0,
"Periodic composites " +
1310 id1s +
" and " + id2s +
1314 for (i = 0; i < compOrder[id1].size(); ++i)
1316 int eId1 = compOrder[id1][i];
1317 int eId2 = compOrder[id2][i];
1319 ASSERTL0(compPairs.count(eId1) == 0,
"Already paired.");
1321 if (compPairs.count(eId2) != 0)
1323 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1325 compPairs[eId1] = eId2;
1331 for (i = 0; i < 2; ++i)
1338 if (c[i]->m_geomVec.size() > 0)
1340 for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1342 locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1348 for (
auto &pIt : compPairs)
1350 int ids[2] = {pIt.first, pIt.second};
1351 bool local[2] = {locEdges.count(pIt.first) > 0,
1352 locEdges.count(pIt.second) > 0};
1354 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1355 orientMap.count(ids[1]) > 0,
1356 "Unable to find edge in orientation map.");
1358 allCompPairs[pIt.first] = pIt.second;
1359 allCompPairs[pIt.second] = pIt.first;
1361 for (i = 0; i < 2; ++i)
1368 int other = (i + 1) % 2;
1371 orientMap[ids[i]] == orientMap[ids[other]]
1379 for (i = 0; i < 2; ++i)
1381 int other = (i + 1) % 2;
1384 orientMap[ids[i]] == orientMap[ids[other]]
1389 vector<int> perVerts1 = vertMap[ids[i]];
1390 vector<int> perVerts2 = vertMap[ids[other]];
1392 map<int, pair<int, bool>> tmpMap;
1395 tmpMap[perVerts1[0]] = make_pair(
1396 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1397 tmpMap[perVerts1[1]] = make_pair(
1398 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1402 tmpMap[perVerts1[0]] = make_pair(
1403 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1404 tmpMap[perVerts1[1]] = make_pair(
1405 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1408 for (
auto &mIt : tmpMap)
1414 auto perIt = periodicVerts.find(mIt.first);
1416 if (perIt == periodicVerts.end())
1418 periodicVerts[mIt.first].push_back(ent2);
1419 perIt = periodicVerts.find(mIt.first);
1424 for (j = 0; j < perIt->second.size(); ++j)
1426 if (perIt->second[j].id == mIt.second.first)
1435 perIt->second.push_back(ent2);
1446 for (cnt = i = 0; i < totEdges; ++i)
1448 int edgeId = edgeIds[i];
1450 ASSERTL0(allCompPairs.count(edgeId) > 0,
1451 "Unable to find matching periodic edge.");
1453 int perEdgeId = allCompPairs[edgeId];
1455 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1457 int vId = traceIds[cnt];
1459 auto perId = periodicVerts.find(vId);
1461 if (perId == periodicVerts.end())
1470 orientMap[edgeId] == orientMap[perEdgeId]
1471 ? vertMap[perEdgeId][(j + 1) % 2]
1472 : vertMap[perEdgeId][j];
1476 locVerts.count(perVertexId) > 0);
1478 periodicVerts[vId].push_back(ent);
1485 for (
auto &perIt : periodicVerts)
1488 for (i = 0; i < perIt.second.size(); ++i)
1490 auto perIt2 = periodicVerts.find(perIt.second[i].id);
1491 ASSERTL0(perIt2 != periodicVerts.end(),
1492 "Couldn't find periodic vertex.");
1494 for (j = 0; j < perIt2->second.size(); ++j)
1496 if (perIt2->second[j].id == perIt.first)
1502 for (k = 0; k < perIt.second.size(); ++k)
1504 if (perIt2->second[j].id == perIt.second[k].id)
1513 perIt.second.push_back(perIt2->second[j]);
1521 for (
auto &perIt : periodicVerts)
1523 if (locVerts.count(perIt.first) > 0)
1533 m_graph->GetCompositeOrdering();
1535 m_graph->GetBndRegionOrdering();
1550 map<int, RotPeriodicInfo> rotComp;
1551 map<int, int> perComps;
1552 map<int, vector<int>> allVerts;
1553 map<int, SpatialDomains::PointGeomVector> allCoord;
1554 map<int, vector<int>> allEdges;
1555 map<int, vector<StdRegions::Orientation>> allOrient;
1560 int region1ID, region2ID, i, j, k, cnt;
1564 for (i = 0; i < (*m_exp).size(); ++i)
1566 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1568 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1569 locVerts.insert(
id);
1572 for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1574 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1575 locEdges.insert(
id);
1582 for (
auto &it : bregions)
1588 if (locBCond->GetBoundaryConditionType() !=
1595 region1ID = it.first;
1597 std::static_pointer_cast<
1599 ->m_connectedBoundaryRegion;
1603 "Boundary region " +
1604 boost::lexical_cast<string>(region1ID) +
1605 " should only contain 1 composite.");
1613 if (vComm->GetSize() == 1)
1615 cId1 = it.second->begin()->first;
1616 cId2 = bregions.find(region2ID)->second->begin()->first;
1620 cId1 = bndRegOrder.find(region1ID)->second[0];
1621 cId2 = bndRegOrder.find(region2ID)->second[0];
1625 if (boost::icontains(locBCond->GetUserDefined(),
"Rotated"))
1627 vector<string> tmpstr;
1629 boost::split(tmpstr, locBCond->GetUserDefined(),
1630 boost::is_any_of(
":"));
1632 if (boost::iequals(tmpstr[0],
"Rotated"))
1635 "Expected Rotated user defined string to "
1636 "contain direction and rotation angle "
1637 "and optionally a tolerance, "
1638 "i.e. Rotated:dir:PI/2:1e-6");
1640 ASSERTL1((tmpstr[1] ==
"x") || (tmpstr[1] ==
"y") ||
1643 "not specified as x,y or z");
1646 RotInfo.
m_dir = (tmpstr[1] ==
"x") ? 0
1647 : (tmpstr[1] ==
"y") ? 1
1654 if (tmpstr.size() == 4)
1659 boost::lexical_cast<NekDouble>(tmpstr[3]);
1664 "failed to cast tolerance input "
1665 "to a double value in Rotated"
1666 "boundary information");
1671 RotInfo.
m_tol = 1e-8;
1673 rotComp[cId1] = RotInfo;
1678 it.second->begin()->second;
1680 vector<unsigned int> tmpOrder;
1688 for (i = 0; i < c->m_geomVec.size(); ++i)
1691 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1693 ASSERTL1(faceGeom,
"Unable to cast to shared ptr");
1696 int faceId = c->m_geomVec[i]->GetGlobalID();
1697 locFaces.insert(faceId);
1701 if (vComm->GetSize() == 1)
1703 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1708 vector<int> vertList, edgeList;
1710 vector<StdRegions::Orientation> orientVec;
1711 for (j = 0; j < faceGeom->GetNumVerts(); ++j)
1713 vertList.push_back(faceGeom->GetVid(j));
1714 edgeList.push_back(faceGeom->GetEid(j));
1715 coordVec.push_back(faceGeom->GetVertex(j));
1716 orientVec.push_back(faceGeom->GetEorient(j));
1719 allVerts[faceId] = vertList;
1720 allEdges[faceId] = edgeList;
1721 allCoord[faceId] = coordVec;
1722 allOrient[faceId] = orientVec;
1727 if (vComm->GetSize() == 1)
1729 compOrder[it.second->begin()->first] = tmpOrder;
1735 if (perComps.count(cId1) == 0)
1737 if (perComps.count(cId2) == 0)
1739 perComps[cId1] = cId2;
1743 std::stringstream ss;
1744 ss <<
"Boundary region " << cId2 <<
" should be "
1745 <<
"periodic with " << perComps[cId2] <<
" but "
1746 <<
"found " << cId1 <<
" instead!";
1747 ASSERTL0(perComps[cId2] == cId1, ss.str());
1752 std::stringstream ss;
1753 ss <<
"Boundary region " << cId1 <<
" should be "
1754 <<
"periodic with " << perComps[cId1] <<
" but "
1755 <<
"found " << cId2 <<
" instead!";
1756 ASSERTL0(perComps[cId1] == cId1, ss.str());
1763 for (
auto &cIt : perComps)
1765 idmax = max(idmax, cIt.first);
1766 idmax = max(idmax, cIt.second);
1771 for (
auto &cIt : perComps)
1773 perid[cIt.first] = cIt.second;
1777 for (
int i = 0; i < idmax; ++i)
1783 if (perComps.count(perid[i]))
1787 perComps[i] = perid[i];
1794 int n = vComm->GetSize();
1795 int p = vComm->GetRank();
1805 rotcounts[
p] = rotComp.size();
1811 for (i = 1; i < n; ++i)
1813 rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
1822 auto rIt = rotComp.begin();
1824 for (i = 0; rIt != rotComp.end(); ++rIt)
1826 compid[rotoffset[
p] + i] = rIt->first;
1827 rotdir[rotoffset[
p] + i] = rIt->second.m_dir;
1828 rotangle[rotoffset[
p] + i] = rIt->second.m_angle;
1829 rottol[rotoffset[
p] + i++] = rIt->second.m_tol;
1838 for (i = 0; i < totrot; ++i)
1842 rotComp[compid[i]] = rinfo;
1847 facecounts[
p] = locFaces.size();
1853 for (i = 1; i < n; ++i)
1855 faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
1869 auto sIt = locFaces.begin();
1870 for (i = 0; sIt != locFaces.end(); ++sIt)
1872 faceIds[faceoffset[
p] + i] = *sIt;
1873 faceVerts[faceoffset[
p] + i++] = allVerts[*sIt].size();
1896 for (i = 0; i < n; ++i)
1898 if (facecounts[i] > 0)
1901 faceVerts + faceoffset[i], 1);
1912 for (i = 1; i < n; ++i)
1914 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1928 for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
1930 for (j = 0; j < allVerts[*sIt].size(); ++j)
1932 int vertId = allVerts[*sIt][j];
1933 vertIds[vertoffset[
p] + cnt] = vertId;
1934 vertX[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(0);
1935 vertY[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(1);
1936 vertZ[vertoffset[
p] + cnt] = (*allCoord[*sIt][j])(2);
1937 edgeIds[vertoffset[
p] + cnt] = allEdges[*sIt][j];
1938 edgeOrt[vertoffset[
p] + cnt++] = allOrient[*sIt][j];
1953 map<int, vector<int>> vertMap;
1954 map<int, vector<int>> edgeMap;
1955 map<int, SpatialDomains::PointGeomVector> coordMap;
1961 map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
1962 map<int, pair<int, int>> eIdMap;
1964 for (cnt = i = 0; i < totFaces; ++i)
1966 vector<int> edges(faceVerts[i]);
1967 vector<int> verts(faceVerts[i]);
1973 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1975 edges[j] = edgeIds[cnt];
1976 verts[j] = vertIds[cnt];
1980 vCoMap[vertIds[cnt]] = coord[j];
1983 auto testIns = eIdMap.insert(make_pair(
1985 make_pair(vertIds[tmp + j],
1986 vertIds[tmp + ((j + 1) % faceVerts[i])])));
1988 if (testIns.second ==
false)
2011 swap(testIns.first->second.first,
2012 testIns.first->second.second);
2016 vertMap[faceIds[i]] = verts;
2017 edgeMap[faceIds[i]] = edges;
2018 coordMap[faceIds[i]] = coord;
2036 map<int, map<StdRegions::Orientation, vector<int>>> vmap;
2037 map<int, map<StdRegions::Orientation, vector<int>>> emap;
2039 map<StdRegions::Orientation, vector<int>> quadVertMap;
2049 map<StdRegions::Orientation, vector<int>> quadEdgeMap;
2059 map<StdRegions::Orientation, vector<int>> triVertMap;
2063 map<StdRegions::Orientation, vector<int>> triEdgeMap;
2067 vmap[3] = triVertMap;
2068 vmap[4] = quadVertMap;
2069 emap[3] = triEdgeMap;
2070 emap[4] = quadEdgeMap;
2072 map<int, int> allCompPairs;
2076 map<int, int> fIdToCompId;
2081 for (
auto &cIt : perComps)
2084 const int id1 = cIt.first;
2085 const int id2 = cIt.second;
2086 std::string id1s = boost::lexical_cast<string>(id1);
2087 std::string id2s = boost::lexical_cast<string>(id2);
2089 if (compMap.count(id1) > 0)
2091 c[0] = compMap[id1];
2094 if (compMap.count(id2) > 0)
2096 c[1] = compMap[id2];
2102 map<int, int> compPairs;
2105 "Unable to find composite " + id1s +
" in order map.");
2107 "Unable to find composite " + id2s +
" in order map.");
2108 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
2109 "Periodic composites " + id1s +
" and " + id2s +
2110 " should have the same number of elements.");
2111 ASSERTL0(compOrder[id1].size() > 0,
"Periodic composites " +
2112 id1s +
" and " + id2s +
2116 for (i = 0; i < compOrder[id1].size(); ++i)
2118 int eId1 = compOrder[id1][i];
2119 int eId2 = compOrder[id2][i];
2121 ASSERTL0(compPairs.count(eId1) == 0,
"Already paired.");
2124 if (compPairs.count(eId2) != 0)
2126 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
2128 compPairs[eId1] = eId2;
2131 fIdToCompId[eId1] = id1;
2132 fIdToCompId[eId2] = id2;
2138 for (
auto &pIt : compPairs)
2140 int ids[2] = {pIt.first, pIt.second};
2141 bool local[2] = {locFaces.count(pIt.first) > 0,
2142 locFaces.count(pIt.second) > 0};
2144 ASSERTL0(coordMap.count(ids[0]) > 0 &&
2145 coordMap.count(ids[1]) > 0,
2146 "Unable to find face in coordinate map");
2148 allCompPairs[pIt.first] = pIt.second;
2149 allCompPairs[pIt.second] = pIt.first;
2154 coordMap[ids[0]], coordMap[ids[1]]};
2156 ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
2157 "Two periodic faces have different number "
2166 bool rotbnd =
false;
2173 if (rotComp.count(fIdToCompId[pIt.first]))
2176 dir = rotComp[fIdToCompId[pIt.first]].m_dir;
2177 angle = rotComp[fIdToCompId[pIt.first]].m_angle;
2178 tol = rotComp[fIdToCompId[pIt.first]].m_tol;
2182 for (i = 0; i < 2; ++i)
2190 int other = (i + 1) % 2;
2193 sign = (i == 0) ? 1.0 : -1.0;
2196 if (tmpVec[0].size() == 3)
2199 tmpVec[i], tmpVec[other], rotbnd, dir,
2205 tmpVec[i], tmpVec[other], rotbnd, dir,
2215 int nFaceVerts = vertMap[ids[0]].size();
2218 for (i = 0; i < 2; ++i)
2220 int other = (i + 1) % 2;
2223 sign = (i == 0) ? 1.0 : -1.0;
2226 if (tmpVec[0].size() == 3)
2229 tmpVec[i], tmpVec[other], rotbnd, dir,
2235 tmpVec[i], tmpVec[other], rotbnd, dir,
2239 if (nFaceVerts == 3)
2244 "Unsupported face orientation for face " +
2245 boost::lexical_cast<string>(ids[i]));
2249 vector<int> per1 = vertMap[ids[i]];
2250 vector<int> per2 = vertMap[ids[other]];
2254 map<int, pair<int, bool>> tmpMap;
2258 for (j = 0; j < nFaceVerts; ++j)
2260 int v = vmap[nFaceVerts][o][j];
2262 make_pair(per2[v], locVerts.count(per2[v]) > 0);
2266 for (
auto &mIt : tmpMap)
2273 auto perIt = periodicVerts.find(mIt.first);
2275 if (perIt == periodicVerts.end())
2279 periodicVerts[mIt.first].push_back(ent2);
2280 perIt = periodicVerts.find(mIt.first);
2287 for (k = 0; k < perIt->second.size(); ++k)
2289 if (perIt->second[k].id == mIt.second.first)
2295 if (k == perIt->second.size())
2297 perIt->second.push_back(ent2);
2305 for (i = 0; i < 2; ++i)
2307 int other = (i + 1) % 2;
2310 sign = (i == 0) ? 1.0 : -1.0;
2312 if (tmpVec[0].size() == 3)
2315 tmpVec[i], tmpVec[other], rotbnd, dir,
2321 tmpVec[i], tmpVec[other], rotbnd, dir,
2325 vector<int> per1 = edgeMap[ids[i]];
2326 vector<int> per2 = edgeMap[ids[other]];
2328 map<int, pair<int, bool>> tmpMap;
2330 for (j = 0; j < nFaceVerts; ++j)
2332 int e = emap[nFaceVerts][o][j];
2334 make_pair(per2[e], locEdges.count(per2[e]) > 0);
2337 for (
auto &mIt : tmpMap)
2344 auto perIt = periodicEdges.find(mIt.first);
2346 if (perIt == periodicEdges.end())
2348 periodicEdges[mIt.first].push_back(ent2);
2349 perIt = periodicEdges.find(mIt.first);
2353 for (k = 0; k < perIt->second.size(); ++k)
2355 if (perIt->second[k].id == mIt.second.first)
2361 if (k == perIt->second.size())
2363 perIt->second.push_back(ent2);
2373 ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
2374 "At this point the size of allCompPairs "
2375 "should have been the same as fIdToCompId");
2379 map<int, int> eIdToCompId;
2385 for (cnt = i = 0; i < totFaces; ++i)
2387 bool rotbnd =
false;
2392 int faceId = faceIds[i];
2394 ASSERTL0(allCompPairs.count(faceId) > 0,
2395 "Unable to find matching periodic face. faceId = " +
2396 boost::lexical_cast<string>(faceId));
2398 int perFaceId = allCompPairs[faceId];
2401 ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
2402 "Face " + boost::lexical_cast<string>(faceId) +
2403 " not found in fIdtoCompId map");
2404 if (rotComp.count(fIdToCompId[faceId]))
2407 dir = rotComp[fIdToCompId[faceId]].m_dir;
2408 angle = rotComp[fIdToCompId[faceId]].m_angle;
2409 tol = rotComp[fIdToCompId[faceId]].m_tol;
2412 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
2414 int vId = vertIds[cnt];
2416 auto perId = periodicVerts.find(vId);
2418 if (perId == periodicVerts.end())
2427 coordMap[faceId], coordMap[perFaceId]};
2429 int nFaceVerts = tmpVec[0].size();
2433 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2436 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2442 vertMap[perFaceId][vmap[nFaceVerts][o][j]];
2446 locVerts.count(perVertexId) > 0);
2448 periodicVerts[vId].push_back(ent);
2451 int eId = edgeIds[cnt];
2453 perId = periodicEdges.find(eId);
2459 eIdToCompId[eId] = fIdToCompId[faceId];
2462 if (perId == periodicEdges.end())
2470 coordMap[faceId], coordMap[perFaceId]};
2472 int nFaceEdges = tmpVec[0].size();
2476 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2479 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2485 edgeMap[perFaceId][emap[nFaceEdges][o][j]];
2488 locEdges.count(perEdgeId) > 0);
2490 periodicEdges[eId].push_back(ent);
2496 eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
2504 for (
auto &perIt : periodicVerts)
2507 for (i = 0; i < perIt.second.size(); ++i)
2510 auto perIt2 = periodicVerts.find(perIt.second[i].id);
2511 ASSERTL0(perIt2 != periodicVerts.end(),
2512 "Couldn't find periodic vertex.");
2517 for (j = 0; j < perIt2->second.size(); ++j)
2519 if (perIt2->second[j].id == perIt.first)
2524 for (k = 0; k < perIt.second.size(); ++k)
2526 if (perIt2->second[j].id == perIt.second[k].id)
2532 if (k == perIt.second.size())
2534 perIt.second.push_back(perIt2->second[j]);
2540 for (
auto &perIt : periodicEdges)
2542 for (i = 0; i < perIt.second.size(); ++i)
2544 auto perIt2 = periodicEdges.find(perIt.second[i].id);
2545 ASSERTL0(perIt2 != periodicEdges.end(),
2546 "Couldn't find periodic edge.");
2548 for (j = 0; j < perIt2->second.size(); ++j)
2550 if (perIt2->second[j].id == perIt.first)
2555 for (k = 0; k < perIt.second.size(); ++k)
2557 if (perIt2->second[j].id == perIt.second[k].id)
2563 if (k == perIt.second.size())
2565 perIt.second.push_back(perIt2->second[j]);
2572 for (
auto &perIt : periodicEdges)
2574 bool rotbnd =
false;
2580 auto eIt = eIdMap.find(perIt.first);
2582 *vCoMap[eIt->second.second]};
2585 if (rotComp.count(eIdToCompId[perIt.first]))
2588 dir = rotComp[eIdToCompId[perIt.first]].m_dir;
2589 angle = rotComp[eIdToCompId[perIt.first]].m_angle;
2590 tol = rotComp[eIdToCompId[perIt.first]].m_tol;
2596 for (i = 0; i < perIt.second.size(); ++i)
2598 eIt = eIdMap.find(perIt.second[i].id);
2601 *vCoMap[eIt->second.first],
2602 *vCoMap[eIt->second.second]};
2604 int vMap[2] = {-1, -1};
2610 r.
Rotate(v[0], dir, angle);
2612 if (r.
dist(
w[0]) < tol)
2618 r.
Rotate(v[1], dir, angle);
2619 if (r.
dist(
w[0]) < tol)
2626 "Unable to align rotationally "
2627 "periodic edge vertex");
2634 0.5 * (
w[0](0) - v[0](0) +
w[1](0) - v[1](0));
2636 0.5 * (
w[0](1) - v[0](1) +
w[1](1) - v[1](1));
2638 0.5 * (
w[0](2) - v[0](2) +
w[1](2) - v[1](2));
2640 for (j = 0; j < 2; ++j)
2645 for (k = 0; k < 2; ++k)
2651 if (
sqrt((x1 - x) * (x1 - x) +
2652 (y1 - y) * (y1 - y) +
2653 (z1 -
z) * (z1 -
z)) < 1e-8)
2662 ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
2663 "Unable to align periodic edge vertex.");
2664 ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
2665 (vMap[1] == 0 || vMap[1] == 1) &&
2666 (vMap[0] != vMap[1]),
2667 "Unable to align periodic edge vertex.");
2681 for (
auto &perIt : periodicVerts)
2683 if (locVerts.count(perIt.first) > 0)
2689 for (
auto &perIt : periodicEdges)
2691 if (locEdges.count(perIt.first) > 0)
2699 ASSERTL1(
false,
"Not setup for this expansion");
2711 "Routine currently only tested for HybridDGHelmholtz");
2714 "Full matrix global systems are not supported for HDG "
2718 "The local to global map is not set up for the requested "
2727 (*m_globalBndMat)[mkey] = glo_matrix;
2731 glo_matrix = matrixIter->second;
2746 bool returnval =
true;
2765 int vSame = (returnval ? 1 : 0);
2768 return (vSame == 1);
2783 for (
int v = 0; v < 2; ++v)
2788 if (vertExp->GetLeftAdjacentElementExp()
2791 (*m_exp)[i]->GetGeom()->GetGlobalID())
2918 "Routine not set up for Gauss Lagrange points distribution");
2935 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2947 bool PutFwdInBwdOnBCs,
bool DoExchange)
2979 for (
int n = 0, cnt = 0; n < nexp; ++n)
2984 for (
int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2987 m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2992 exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
3008 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
3019 bool PutFwdInBwdOnBCs)
3032 bool PutFwdInBwdOnBCs)
3036 if (PutFwdInBwdOnBCs)
3039 for (
int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3041 if (bndConditions[n]->GetBoundaryConditionType() ==
3044 auto ne = bndCondExpansions[n]->GetExpSize();
3045 for (
int e = 0; e < ne; ++e)
3047 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3048 auto id2 =
m_trace->GetPhys_Offset(
3049 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3055 else if (bndConditions[n]->GetBoundaryConditionType() ==
3057 bndConditions[n]->GetBoundaryConditionType() ==
3060 auto ne = bndCondExpansions[n]->GetExpSize();
3061 for (
int e = 0; e < ne; ++e)
3063 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3064 auto id2 =
m_trace->GetPhys_Offset(
3065 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3070 else if (bndConditions[n]->GetBoundaryConditionType() !=
3073 ASSERTL0(
false,
"Method not set up for this "
3074 "boundary condition.");
3080 for (
int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3082 if (bndConditions[n]->GetBoundaryConditionType() ==
3085 auto ne = bndCondExpansions[n]->GetExpSize();
3086 for (
int e = 0; e < ne; ++e)
3088 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3089 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3090 auto id2 =
m_trace->GetPhys_Offset(
3091 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3098 else if (bndConditions[n]->GetBoundaryConditionType() ==
3100 bndConditions[n]->GetBoundaryConditionType() ==
3104 for (
int e = 0; e < ne; ++e)
3106 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3107 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3110 "Method not set up for non-zero "
3111 "Neumann boundary condition");
3112 auto id2 =
m_trace->GetPhys_Offset(
3113 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3120 else if (bndConditions[n]->GetBoundaryConditionType() ==
3125 else if (bndConditions[n]->GetBoundaryConditionType() !=
3129 "Method not set up for this boundary "
3169 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3204 m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3210 int n,
p, offset, phys_offset;
3215 "input array is of insufficient length");
3216 for (n = 0; n < nexp; ++n)
3219 for (
p = 0;
p < (*m_exp)[n]->GetNtraces(); ++
p)
3222 m_trace->GetPhys_Offset(elmtToTrace[n][
p]->GetElmtId());
3223 (*m_exp)[n]->GetTracePhysVals(
p, elmtToTrace[n][
p],
3224 inarray + phys_offset,
3225 t_tmp = outarray + offset);
3253 int n, offset, t_offset;
3260 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3263 if ((*
m_exp)[n]->GetBasis(0)->GetBasisType() !=
3267 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3268 if (negatedFluxNormal[2 * n])
3270 outarray[offset] -= Fn[t_offset];
3274 outarray[offset] += Fn[t_offset];
3277 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3278 if (negatedFluxNormal[2 * n + 1])
3280 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3285 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3293 static int sav_ncoeffs = 0;
3294 if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3304 m_Ixm = BASE->GetI(coords);
3306 m_Ixp = BASE->GetI(coords);
3307 sav_ncoeffs = e_ncoeffs;
3310 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3311 if (negatedFluxNormal[2 * n])
3313 for (j = 0; j < e_ncoeffs; j++)
3315 outarray[offset + j] -=
3316 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3321 for (j = 0; j < e_ncoeffs; j++)
3323 outarray[offset + j] +=
3324 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3328 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3329 if (negatedFluxNormal[2 * n + 1])
3331 for (j = 0; j < e_ncoeffs; j++)
3333 outarray[offset + j] -=
3334 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3339 for (j = 0; j < e_ncoeffs; j++)
3341 outarray[offset + j] +=
3342 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3356 m_trace->IProductWRTBase(Fn, Fcoeffs);
3362 int e, n, offset, t_offset;
3371 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3373 t_offset =
GetTrace()->GetPhys_Offset(
3374 elmtToTrace[n][e]->GetElmtId());
3375 (*m_exp)[n]->AddEdgeNormBoundaryInt(
3376 e, elmtToTrace[n][e], Fn + t_offset,
3377 e_outarray = outarray + offset);
3386 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3388 t_offset =
GetTrace()->GetPhys_Offset(
3389 elmtToTrace[n][e]->GetElmtId());
3390 (*m_exp)[n]->AddFaceNormBoundaryInt(
3391 e, elmtToTrace[n][e], Fn + t_offset,
3392 e_outarray = outarray + offset);
3424 "tested for 1D expansion");
3426 m_trace->IProductWRTBase(Fwd, Coeffs);
3428 m_trace->IProductWRTBase(Bwd, Coeffs);
3438 const bool PhysSpaceForcing)
3440 int i, n, cnt, nbndry;
3448 if (PhysSpaceForcing)
3460 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
3461 int NumDirBCs =
m_traceMap->GetNumLocalDirBndCoeffs();
3468 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
3477 for (cnt = n = 0; n < nexp; ++n)
3479 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3480 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3486 Floc =
Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3490 m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3506 locid = bndCondMap[cnt + j];
3507 loclambda[locid] = Sign[locid] * bndcoeffs[j];
3518 locid = bndCondMap[cnt + j];
3519 bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3536 if (GloBndDofs - NumDirBCs > 0)
3541 LinSys->Solve(bndrhs, loclambda,
m_traceMap);
3559 out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3601 const std::string varName,
3617 npoints = locExpList->GetNpoints();
3622 locExpList->GetCoords(x0, x1, x2);
3646 0, (std::static_pointer_cast<
3649 ->m_dirichletCondition)
3650 .Evaluate(x0[0], x1[0], x2[0], time));
3658 0, (std::static_pointer_cast<
3661 ->m_neumannCondition)
3662 .Evaluate(x0[0], x1[0], x2[0], time));
3668 0, (std::static_pointer_cast<
3672 .Evaluate(x0[0], x1[0], x2[0], time));
3686 "This type of BC not implemented yet");
3695 std::static_pointer_cast<
3700 valuesExp(npoints, 1.0);
3702 string bcfilename = bcPtr->m_filename;
3703 string exprbcs = bcPtr->m_expr;
3705 if (bcfilename !=
"")
3707 locExpList->ExtractCoeffsFromFile(
3708 bcfilename, bcPtr->GetComm(), varName,
3709 locExpList->UpdateCoeffs());
3710 locExpList->BwdTrans(locExpList->GetCoeffs(),
3711 locExpList->UpdatePhys());
3712 valuesFile = locExpList->GetPhys();
3718 std::static_pointer_cast<
3721 ->m_dirichletCondition;
3723 condition.
Evaluate(x0, x1, x2, time, valuesExp);
3727 locExpList->UpdatePhys(), 1);
3728 locExpList->FwdTransBndConstrained(
3729 locExpList->GetPhys(), locExpList->UpdateCoeffs());
3735 std::static_pointer_cast<
3740 if (bcfilename !=
"")
3742 locExpList->ExtractCoeffsFromFile(
3743 bcfilename, bcPtr->GetComm(), varName,
3744 locExpList->UpdateCoeffs());
3745 locExpList->BwdTrans(locExpList->GetCoeffs(),
3746 locExpList->UpdatePhys());
3751 std::static_pointer_cast<
3754 ->m_neumannCondition;
3755 condition.
Evaluate(x0, x1, x2, time,
3756 locExpList->UpdatePhys());
3759 locExpList->IProductWRTBase(locExpList->GetPhys(),
3760 locExpList->UpdateCoeffs());
3766 std::static_pointer_cast<
3772 if (bcfilename !=
"")
3774 locExpList->ExtractCoeffsFromFile(
3775 bcfilename, bcPtr->GetComm(), varName,
3776 locExpList->UpdateCoeffs());
3777 locExpList->BwdTrans(locExpList->GetCoeffs(),
3778 locExpList->UpdatePhys());
3783 std::static_pointer_cast<
3787 condition.
Evaluate(x0, x1, x2, time,
3788 locExpList->UpdatePhys());
3791 locExpList->IProductWRTBase(locExpList->GetPhys(),
3792 locExpList->UpdateCoeffs());
3802 "This type of BC not implemented yet");
3831 id2 =
m_trace->GetPhys_Offset(
3832 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3847 id2 =
m_trace->GetPhys_Offset(
3848 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3859 "Method not set up for this boundary condition.");
3876 map<int, int> VertGID;
3901 VertGID[Vid] = cnt++;
3910 for (i = 0; i < exp->GetNverts(); ++i)
3912 id = exp->GetGeom()->GetVid(i);
3914 if (VertGID.count(
id) > 0)
3916 bid = VertGID.find(
id)->second;
3918 "Edge already set");
3926 "Failed to visit all boundary condtiions");
3931 map<int, int> globalIdMap;
3940 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3966 m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3969 globalIdMap[(*tmp)[0].first->GetGlobalID()];
3977 map<int, int> globalIdMap;
3987 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
4011 m_graph->GetElementsFromFace(exp2d->GetGeom2D());
4013 globalIdMap[tmp->at(0).first->GetGlobalID()];
4020 ASSERTL1(
false,
"Not setup for this expansion");
4030 std::shared_ptr<ExpList> &result,
4031 const bool DeclareCoeffPhysArrays)
4034 std::vector<unsigned int> eIDs;
4040 for (cnt = n = 0; n < i; ++n)
4048 eIDs.push_back(ElmtID[cnt + n]);
4055 if (DeclareCoeffPhysArrays)
4058 int offsetOld, offsetNew;
4060 for (n = 0; n < result->GetExpSize(); ++n)
4062 nq =
GetExp(ElmtID[cnt + n])->GetTotPoints();
4064 offsetNew = result->GetPhys_Offset(n);
4066 tmp2 = result->UpdatePhys() + offsetNew, 1);
4067 nq =
GetExp(ElmtID[cnt + n])->GetNcoeffs();
4069 offsetNew = result->GetCoeff_Offset(n);
4071 tmp2 = result->UpdateCoeffs() + offsetNew, 1);
4107 map<int, RobinBCInfoSharedPtr> returnval;
4123 int npoints = locExpList->GetNpoints();
4129 locExpList->GetCoords(x0, x1, x2);
4132 std::static_pointer_cast<
4134 ->m_robinPrimitiveCoeff;
4137 coeffeqn.
Evaluate(x0, x1, x2, 0.0, coeffphys);
4139 for (e = 0; e < locExpList->GetExpSize(); ++e)
4144 Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4146 elmtid = ElmtID[cnt + e];
4148 if (returnval.count(elmtid) != 0)
4150 rInfo->next = returnval.find(elmtid)->second;
4152 returnval[elmtid] = rInfo;
4175 int i, e, ncoeff_edge;
4185 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4189 edge_lambda = loc_lambda;
4197 int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4200 for (e = 0; e < nEdges; ++e)
4202 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4203 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4205 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4206 elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4208 edge_lambda = edge_lambda + ncoeff_edge;
4211 (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs +
m_coeff_offset[i],
4212 elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4213 cnt += (*m_exp)[i]->GetNcoeffs();
4246 int i, cnt, e, ncoeff_trace;
4253 int nq_elmt, nm_elmt;
4254 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4259 trace_lambda = loc_lambda;
4263 int num_points[] = {0, 0, 0};
4264 int num_modes[] = {0, 0, 0};
4269 nq_elmt = (*m_exp)[i]->GetTotPoints();
4270 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4274 out_tmp = force + nm_elmt;
4277 for (
int j = 0; j < dim; ++j)
4279 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4280 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4286 int nTraces = (*m_exp)[i]->GetNtraces();
4289 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4291 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4292 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4294 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4297 elmtToTrace[i][e]->SetCoeffsToOrientation(
4298 edgedir, traceCoeffs[e], traceCoeffs[e]);
4304 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4306 trace_lambda = trace_lambda + ncoeff_trace;
4321 num_modes[0], PkeyQ1);
4323 num_modes[1], PkeyQ2);
4325 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4326 (*
m_exp)[i]->GetGeom());
4328 BkeyQ1, BkeyQ2, qGeom);
4336 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4338 num_modes[0], PkeyT1);
4340 num_modes[1], PkeyT2);
4342 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4343 (*
m_exp)[i]->GetGeom());
4345 BkeyT1, BkeyT2, tGeom);
4357 num_modes[0], PkeyH1);
4359 num_modes[1], PkeyH2);
4361 num_modes[2], PkeyH3);
4363 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4364 (*
m_exp)[i]->GetGeom());
4366 BkeyH1, BkeyH2, BkeyH3, hGeom);
4374 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4376 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4378 num_modes[0], PkeyT1);
4380 num_modes[1], PkeyT2);
4382 num_modes[2], PkeyT3);
4384 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4385 (*
m_exp)[i]->GetGeom());
4387 BkeyT1, BkeyT2, BkeyT3, tGeom);
4392 "Wrong shape type, HDG postprocessing is not "
4399 elmtToTrace[i], traceCoeffs, out_tmp);
4400 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4401 ppExp->IProductWRTDerivBase(0, qrhs, force);
4405 elmtToTrace[i], traceCoeffs, out_tmp);
4407 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4408 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4410 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4413 (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs +
m_coeff_offset[i], qrhs);
4414 force[0] = (*m_exp)[i]->Integral(qrhs);
4419 ppExp->DetShapeType(), *ppExp);
4425 out = (*lapsys) * in;
4429 ppExp->BwdTrans(out.
GetPtr(), work);
4430 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray +
m_coeff_offset[i]);
4459 m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4461 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 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
Exapnsion 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< 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.