52 namespace MultiRegions
64 DisContField::DisContField()
65 :
ExpList(), m_bndCondExpansions(), m_bndConditions(),
84 const std::string &variable,
const bool SetUpJustDG,
85 const bool DeclareCoeffPhysArrays,
87 const std::string bcvariable)
88 :
ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
92 if (bcvariable ==
"NotSet")
101 if (bcvar.compare(
"DefaultVar") != 0)
106 DeclareCoeffPhysArrays);
107 if (DeclareCoeffPhysArrays)
130 for (
int e = 0; e < locExpList->GetExpSize(); ++e)
132 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
133 locExpList->GetExp(e));
134 locExpList->GetExp(e)->SetAdjacentElementExp(
135 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
140 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
142 std::string ProjectStr =
m_session->GetSolverInfo(
"PROJECTION");
143 if ((ProjectStr ==
"MixedCGDG") ||
144 (ProjectStr ==
"Mixed_CG_Discontinuous"))
187 if (
m_session->DefinesCmdLineArgument(
"verbose"))
200 for (
int i = 0; i <
m_exp->size(); ++i)
202 for (
int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
207 exp->SetAdjacentElementExp(j, (*
m_exp)[i]);
208 (*m_exp)[i]->SetTraceExp(j, exp);
226 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
233 std::unordered_map<int, pair<int, int>> perTraceToExpMap;
234 for (cnt = n = 0; n <
m_exp->size(); ++n)
236 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
238 auto it = periodicTraces.find((*
m_exp)[n]->GetGeom()->GetTid(v));
240 if (it != periodicTraces.end())
242 perTraceToExpMap[it->first] = make_pair(n, v);
251 for (cnt = n = 0; n <
m_exp->size(); ++n)
253 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
260 for (cnt = n = 0; n <
m_exp->size(); ++n)
262 for (
int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
264 int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
267 auto it = periodicTraces.find(GeomId);
269 if (it != periodicTraces.end())
272 auto it2 = perTraceToExpMap.find(ent.
id);
274 if (it2 == perTraceToExpMap.end())
276 if (
m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
283 ASSERTL1(
false,
"Periodic trace not found!");
288 "Periodic trace in non-forward space?");
290 int offset =
m_trace->GetPhys_Offset(
291 (
m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
293 int offset2 =
m_trace->GetPhys_Offset(
294 (
m_traceMap->GetElmtToTrace())[it2->second.first]
310 int nquad = elmtToTrace[n][v]->GetNumPoints(0);
312 vector<int> tmpBwd(nquad);
313 vector<int> tmpFwd(nquad);
317 for (
int i = 0; i < nquad; ++i)
319 tmpBwd[i] = offset2 + i;
320 tmpFwd[i] = offset + i;
325 for (
int i = 0; i < nquad; ++i)
327 tmpBwd[i] = offset2 + i;
328 tmpFwd[i] = offset + nquad - i - 1;
332 for (
int i = 0; i < nquad; ++i)
343 int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
344 int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
346 vector<int> tmpBwd(nquad1 * nquad2);
347 vector<int> tmpFwd(nquad1 * nquad2);
357 for (
int i = 0; i < nquad2; ++i)
359 for (
int j = 0; j < nquad1; ++j)
361 tmpBwd[i * nquad1 + j] =
362 offset2 + i * nquad1 + j;
363 tmpFwd[i * nquad1 + j] =
364 offset + j * nquad2 + i;
370 for (
int i = 0; i < nquad2; ++i)
372 for (
int j = 0; j < nquad1; ++j)
374 tmpBwd[i * nquad1 + j] =
375 offset2 + i * nquad1 + j;
376 tmpFwd[i * nquad1 + j] =
377 offset + i * nquad1 + j;
391 for (
int i = 0; i < nquad2; ++i)
393 for (
int j = 0; j < nquad1 / 2; ++j)
395 swap(tmpFwd[i * nquad1 + j],
396 tmpFwd[i * nquad1 + nquad1 - j - 1]);
410 for (
int j = 0; j < nquad1; ++j)
412 for (
int i = 0; i < nquad2 / 2; ++i)
414 swap(tmpFwd[i * nquad1 + j],
415 tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
420 for (
int i = 0; i < nquad1 * nquad2; ++i)
483 if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
484 traceEl->GetRightAdjacentElementTrace() == -1)
497 else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
498 traceEl->GetRightAdjacentElementTrace() != -1)
501 fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*
m_exp)[n].get();
505 ASSERTL2(
false,
"Unconnected trace element!");
517 const std::string &variable)
519 boost::ignore_unused(variable);
525 map<int, int> GeometryToRegionsMap;
533 for (
auto &it : bregions)
535 for (
auto &bregionIt : *it.second)
539 int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
540 GeometryToRegionsMap[id] = it.first;
544 map<int, SpatialDomains::GeometrySharedPtr> EndOfDomain;
547 for (
auto &domIt : domain)
550 for (
int i = 0; i < geomvector->m_geomVec.size(); ++i)
552 for (
int j = 0; j < 2; ++j)
554 int vid = geomvector->m_geomVec[i]->GetVid(j);
555 if (EndOfDomain.count(vid) == 0)
557 EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
561 EndOfDomain.erase(vid);
566 ASSERTL1(EndOfDomain.size() == 2,
"Did not find two ends of domain");
568 for (
auto ®It : EndOfDomain)
570 if (GeometryToRegionsMap.count(regIt.first) != 0)
573 auto iter = GeometryToRegionsMap.find(regIt.first);
574 ASSERTL1(iter != GeometryToRegionsMap.end(),
575 "Failied to find GeometryToRegionMap");
577 int regionId = iter->second;
578 auto bregionsIter = bregions.find(regionId);
579 ASSERTL1(bregionsIter != bregions.end(),
580 "Failed to find boundary region");
583 returnval->AddBoundaryRegions(regionId, breg);
585 auto bconditionsIter = bconditions.find(regionId);
586 ASSERTL1(bconditionsIter != bconditions.end(),
587 "Failed to find boundary collection");
590 bconditionsIter->second;
591 returnval->AddBoundaryConditions(regionId, bcond);
609 const std::string &variable,
611 bool SetToOneSpaceDimension,
613 :
ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
616 if (variable.compare(
"DefaultVar") != 0)
635 const bool DeclareCoeffPhysArrays)
636 :
ExpList(In, DeclareCoeffPhysArrays),
637 m_bndCondExpansions(In.m_bndCondExpansions),
638 m_bndConditions(In.m_bndConditions), m_globalBndMat(In.m_globalBndMat),
639 m_traceMap(In.m_traceMap), m_boundaryTraces(In.m_boundaryTraces),
640 m_periodicVerts(In.m_periodicVerts),
641 m_periodicFwdCopy(In.m_periodicFwdCopy),
642 m_periodicBwdCopy(In.m_periodicBwdCopy),
643 m_leftAdjacentTraces(In.m_leftAdjacentTraces),
644 m_locTraceToTraceMap(In.m_locTraceToTraceMap)
649 *In.
m_trace, DeclareCoeffPhysArrays);
659 const std::string &variable,
const bool SetUpJustDG,
660 const bool DeclareCoeffPhysArrays)
661 :
ExpList(In, DeclareCoeffPhysArrays)
668 if (variable.compare(
"DefaultVar") != 0)
673 if (DeclareCoeffPhysArrays)
701 for (e = 0; e < locExpList->GetExpSize(); ++e)
703 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
704 TraceID[cnt + e], locExpList->GetExp(e));
705 locExpList->GetExp(e)->SetAdjacentElementExp(
706 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
712 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
714 std::string ProjectStr =
717 if ((ProjectStr ==
"MixedCGDG") ||
718 (ProjectStr ==
"Mixed_CG_Discontinuous"))
747 if (SetUpJustDG ==
false)
761 for (e = 0; e < locExpList->GetExpSize(); ++e)
763 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
764 TraceID[cnt + e], locExpList->GetExp(e));
765 locExpList->GetExp(e)->SetAdjacentElementExp(
766 TraceID[cnt + e], (*
m_exp)[ElmtID[cnt + e]]);
812 const bool DeclareCoeffPhysArrays)
830 for (
auto &it : bregions)
835 m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
836 false, bc->GetComm());
847 boost::iequals(type,
"I") || boost::iequals(type,
"CalcBC"))
876 int i, region1ID, region2ID;
879 map<int, int> BregionToVertMap;
883 for (
auto &it : bregions)
888 if (locBCond->GetBoundaryConditionType() !=
894 it.second->begin()->second->m_geomVec[0]->GetGlobalID();
896 BregionToVertMap[it.first] = id;
901 int n = vComm->GetSize();
902 int p = vComm->GetRank();
905 nregions[
p] = BregionToVertMap.size();
912 for (i = 1; i < n; ++i)
914 regOffset[i] = regOffset[i - 1] + nregions[i - 1];
921 for (
auto &iit : BregionToVertMap)
923 bregid[i] = iit.first;
924 bregmap[i++] = iit.second;
925 islocal.insert(iit.first);
931 for (
int i = 0; i < totRegions; ++i)
933 BregionToVertMap[bregid[i]] = bregmap[i];
937 for (
auto &it : bregions)
942 if (locBCond->GetBoundaryConditionType() !=
949 region1ID = it.first;
951 std::static_pointer_cast<
953 ->m_connectedBoundaryRegion;
955 ASSERTL0(BregionToVertMap.count(region1ID) != 0,
956 "Cannot determine vertex of region1ID");
958 ASSERTL0(BregionToVertMap.count(region2ID) != 0,
959 "Cannot determine vertex of region2ID");
963 islocal.count(region2ID) != 0);
971 int region1ID, region2ID, i, j, k, cnt;
975 m_graph->GetCompositeOrdering();
977 m_graph->GetBndRegionOrdering();
984 map<int, int> perComps;
985 map<int, vector<int>> allVerts;
987 map<int, pair<int, StdRegions::Orientation>> allEdges;
990 for (i = 0; i < (*m_exp).size(); ++i)
992 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
994 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1000 for (
auto &it : bregions)
1005 if (locBCond->GetBoundaryConditionType() !=
1012 region1ID = it.first;
1014 std::static_pointer_cast<
1016 ->m_connectedBoundaryRegion;
1021 if (vComm->GetSize() == 1)
1023 cId1 = it.second->begin()->first;
1024 cId2 = bregions.find(region2ID)->second->begin()->first;
1028 cId1 = bndRegOrder.find(region1ID)->second[0];
1029 cId2 = bndRegOrder.find(region2ID)->second[0];
1033 "Boundary region " +
1034 boost::lexical_cast<string>(region1ID) +
1035 " should only contain 1 composite.");
1037 vector<unsigned int> tmpOrder;
1042 it.second->begin()->second;
1044 for (i = 0; i < c->m_geomVec.size(); ++i)
1047 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1049 ASSERTL0(segGeom,
"Unable to cast to shared ptr");
1052 m_graph->GetElementsFromEdge(segGeom);
1054 "The periodic boundaries belong to "
1055 "more than one element of the mesh");
1058 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1061 allEdges[c->m_geomVec[i]->GetGlobalID()] =
1062 make_pair(elmt->at(0).second,
1063 geom->GetEorient(elmt->at(0).second));
1067 if (vComm->GetSize() == 1)
1069 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1072 vector<int> vertList(2);
1073 vertList[0] = segGeom->GetVid(0);
1074 vertList[1] = segGeom->GetVid(1);
1075 allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1078 if (vComm->GetSize() == 1)
1080 compOrder[it.second->begin()->first] = tmpOrder;
1085 if (perComps.count(cId1) == 0)
1087 if (perComps.count(cId2) == 0)
1089 perComps[cId1] = cId2;
1093 std::stringstream ss;
1094 ss <<
"Boundary region " << cId2 <<
" should be "
1095 <<
"periodic with " << perComps[cId2] <<
" but "
1096 <<
"found " << cId1 <<
" instead!";
1097 ASSERTL0(perComps[cId2] == cId1, ss.str());
1102 std::stringstream ss;
1103 ss <<
"Boundary region " << cId1 <<
" should be "
1104 <<
"periodic with " << perComps[cId1] <<
" but "
1105 <<
"found " << cId2 <<
" instead!";
1106 ASSERTL0(perComps[cId1] == cId1, ss.str());
1113 for (
auto &cIt : perComps)
1115 idmax = max(idmax, cIt.first);
1116 idmax = max(idmax, cIt.second);
1121 for (
auto &cIt : perComps)
1123 perid[cIt.first] = cIt.second;
1127 for (
int i = 0; i < idmax; ++i)
1133 if (perComps.count(perid[i]))
1137 perComps[i] = perid[i];
1143 int n = vComm->GetSize();
1144 int p = vComm->GetRank();
1150 edgecounts[
p] = allEdges.size();
1154 for (i = 1; i < n; ++i)
1156 edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
1165 auto sIt = allEdges.begin();
1167 for (i = 0; sIt != allEdges.end(); ++sIt)
1169 edgeIds[edgeoffset[
p] + i] = sIt->first;
1170 edgeIdx[edgeoffset[
p] + i] = sIt->second.first;
1171 edgeOrient[edgeoffset[
p] + i] = sIt->second.second;
1172 edgeVerts[edgeoffset[
p] + i++] = allVerts[sIt->first].size();
1195 for (i = 0; i < n; ++i)
1197 if (edgecounts[i] > 0)
1200 edgeVerts + edgeoffset[i], 1);
1209 for (i = 1; i < n; ++i)
1211 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1215 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1217 for (j = 0; j < allVerts[sIt->first].size(); ++j)
1219 traceIds[vertoffset[
p] + i++] = allVerts[sIt->first][j];
1226 map<int, StdRegions::Orientation> orientMap;
1227 map<int, vector<int>> vertMap;
1229 for (cnt = i = 0; i < totEdges; ++i)
1231 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1232 "Already found edge in orientation map!");
1252 orientMap[edgeIds[i]] = o;
1254 vector<int> verts(edgeVerts[i]);
1256 for (j = 0; j < edgeVerts[i]; ++j)
1258 verts[j] = traceIds[cnt++];
1260 vertMap[edgeIds[i]] = verts;
1267 map<int, int> allCompPairs;
1276 for (
auto &cIt : perComps)
1279 const int id1 = cIt.first;
1280 const int id2 = cIt.second;
1281 std::string id1s = boost::lexical_cast<string>(id1);
1282 std::string id2s = boost::lexical_cast<string>(id2);
1284 if (compMap.count(id1) > 0)
1286 c[0] = compMap[id1];
1289 if (compMap.count(id2) > 0)
1291 c[1] = compMap[id2];
1297 map<int, int> compPairs;
1300 "Unable to find composite " + id1s +
" in order map.");
1302 "Unable to find composite " + id2s +
" in order map.");
1303 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1304 "Periodic composites " + id1s +
" and " + id2s +
1305 " should have the same number of elements.");
1306 ASSERTL0(compOrder[id1].size() > 0,
"Periodic composites " +
1307 id1s +
" and " + id2s +
1311 for (i = 0; i < compOrder[id1].size(); ++i)
1313 int eId1 = compOrder[id1][i];
1314 int eId2 = compOrder[id2][i];
1316 ASSERTL0(compPairs.count(eId1) == 0,
"Already paired.");
1318 if (compPairs.count(eId2) != 0)
1320 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1322 compPairs[eId1] = eId2;
1328 for (i = 0; i < 2; ++i)
1335 if (c[i]->m_geomVec.size() > 0)
1337 for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1339 locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1345 for (
auto &pIt : compPairs)
1347 int ids[2] = {pIt.first, pIt.second};
1348 bool local[2] = {locEdges.count(pIt.first) > 0,
1349 locEdges.count(pIt.second) > 0};
1351 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1352 orientMap.count(ids[1]) > 0,
1353 "Unable to find edge in orientation map.");
1355 allCompPairs[pIt.first] = pIt.second;
1356 allCompPairs[pIt.second] = pIt.first;
1358 for (i = 0; i < 2; ++i)
1365 int other = (i + 1) % 2;
1368 orientMap[ids[i]] == orientMap[ids[other]]
1376 for (i = 0; i < 2; ++i)
1378 int other = (i + 1) % 2;
1381 orientMap[ids[i]] == orientMap[ids[other]]
1386 vector<int> perVerts1 = vertMap[ids[i]];
1387 vector<int> perVerts2 = vertMap[ids[other]];
1389 map<int, pair<int, bool>> tmpMap;
1392 tmpMap[perVerts1[0]] = make_pair(
1393 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1394 tmpMap[perVerts1[1]] = make_pair(
1395 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1399 tmpMap[perVerts1[0]] = make_pair(
1400 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1401 tmpMap[perVerts1[1]] = make_pair(
1402 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1405 for (
auto &mIt : tmpMap)
1411 auto perIt = periodicVerts.find(mIt.first);
1413 if (perIt == periodicVerts.end())
1415 periodicVerts[mIt.first].push_back(ent2);
1416 perIt = periodicVerts.find(mIt.first);
1421 for (j = 0; j < perIt->second.size(); ++j)
1423 if (perIt->second[j].id == mIt.second.first)
1432 perIt->second.push_back(ent2);
1443 for (cnt = i = 0; i < totEdges; ++i)
1445 int edgeId = edgeIds[i];
1447 ASSERTL0(allCompPairs.count(edgeId) > 0,
1448 "Unable to find matching periodic edge.");
1450 int perEdgeId = allCompPairs[edgeId];
1452 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1454 int vId = traceIds[cnt];
1456 auto perId = periodicVerts.find(vId);
1458 if (perId == periodicVerts.end())
1467 orientMap[edgeId] == orientMap[perEdgeId]
1468 ? vertMap[perEdgeId][(j + 1) % 2]
1469 : vertMap[perEdgeId][j];
1473 locVerts.count(perVertexId) > 0);
1475 periodicVerts[vId].push_back(ent);
1482 for (
auto &perIt : periodicVerts)
1485 for (i = 0; i < perIt.second.size(); ++i)
1487 auto perIt2 = periodicVerts.find(perIt.second[i].id);
1488 ASSERTL0(perIt2 != periodicVerts.end(),
1489 "Couldn't find periodic vertex.");
1491 for (j = 0; j < perIt2->second.size(); ++j)
1493 if (perIt2->second[j].id == perIt.first)
1499 for (k = 0; k < perIt.second.size(); ++k)
1501 if (perIt2->second[j].id == perIt.second[k].id)
1510 perIt.second.push_back(perIt2->second[j]);
1518 for (
auto &perIt : periodicVerts)
1520 if (locVerts.count(perIt.first) > 0)
1530 m_graph->GetCompositeOrdering();
1532 m_graph->GetBndRegionOrdering();
1547 map<int, RotPeriodicInfo> rotComp;
1548 map<int, int> perComps;
1549 map<int, vector<int>> allVerts;
1550 map<int, SpatialDomains::PointGeomVector> allCoord;
1551 map<int, vector<int>> allEdges;
1552 map<int, vector<StdRegions::Orientation>> allOrient;
1557 int region1ID, region2ID, i, j, k, cnt;
1561 for (i = 0; i < (*m_exp).size(); ++i)
1563 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1565 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1566 locVerts.insert(
id);
1569 for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1571 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1572 locEdges.insert(
id);
1579 for (
auto &it : bregions)
1585 if (locBCond->GetBoundaryConditionType() !=
1592 region1ID = it.first;
1594 std::static_pointer_cast<
1596 ->m_connectedBoundaryRegion;
1600 "Boundary region " +
1601 boost::lexical_cast<string>(region1ID) +
1602 " should only contain 1 composite.");
1610 if (vComm->GetSize() == 1)
1612 cId1 = it.second->begin()->first;
1613 cId2 = bregions.find(region2ID)->second->begin()->first;
1617 cId1 = bndRegOrder.find(region1ID)->second[0];
1618 cId2 = bndRegOrder.find(region2ID)->second[0];
1622 if (boost::icontains(locBCond->GetUserDefined(),
"Rotated"))
1624 vector<string> tmpstr;
1626 boost::split(tmpstr, locBCond->GetUserDefined(),
1627 boost::is_any_of(
":"));
1629 if (boost::iequals(tmpstr[0],
"Rotated"))
1632 "Expected Rotated user defined string to "
1633 "contain direction and rotation angle "
1634 "and optionally a tolerance, "
1635 "i.e. Rotated:dir:PI/2:1e-6");
1637 ASSERTL1((tmpstr[1] ==
"x") || (tmpstr[1] ==
"y") ||
1640 "not specified as x,y or z");
1643 RotInfo.
m_dir = (tmpstr[1] ==
"x") ? 0
1644 : (tmpstr[1] ==
"y") ? 1
1651 if (tmpstr.size() == 4)
1656 boost::lexical_cast<NekDouble>(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 = boost::lexical_cast<string>(id1);
2084 std::string id2s = boost::lexical_cast<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 boost::lexical_cast<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 boost::lexical_cast<string>(faceId));
2395 int perFaceId = allCompPairs[faceId];
2398 ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
2399 "Face " + boost::lexical_cast<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())
2833 bool PutFwdInBwdOnBCs,
bool DoExchange)
2865 for (
int n = 0, cnt = 0; n < nexp; ++n)
2870 for (
int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2873 m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2878 exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2894 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2900 bool PutFwdInBwdOnBCs)
2903 if (PutFwdInBwdOnBCs)
2912 for (
int e = 0; e < ne; ++e)
2916 auto id2 =
m_trace->GetPhys_Offset(
2917 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2929 for (
int e = 0; e < ne; ++e)
2933 auto id2 =
m_trace->GetPhys_Offset(
2934 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2942 ASSERTL0(
false,
"Method not set up for this "
2943 "boundary condition.");
2955 for (
int e = 0; e < ne; ++e)
2960 auto id2 =
m_trace->GetPhys_Offset(
2961 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2974 for (
int e = 0; e < ne; ++e)
2980 "Method not set up for non-zero "
2981 "Neumann boundary condition");
2982 auto id2 =
m_trace->GetPhys_Offset(
2983 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2998 "Method not set up for this boundary "
3028 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3065 m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3072 int n,
p, offset, phys_offset;
3079 "input array is of insufficient length");
3081 for (n = 0; n < nexp; ++n)
3085 for (
p = 0;
p < (*m_exp)[n]->GetNtraces(); ++
p)
3088 m_trace->GetPhys_Offset(elmtToTrace[n][
p]->GetElmtId());
3089 (*m_exp)[n]->GetTracePhysVals(
p, elmtToTrace[n][
p],
3090 inarray + phys_offset,
3091 t_tmp = outarray + offset);
3120 int n, offset, t_offset;
3130 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3135 if ((*
m_exp)[n]->GetBasis(0)->GetBasisType() !=
3139 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3140 if (negatedFluxNormal[2 * n])
3142 outarray[offset] -= Fn[t_offset];
3146 outarray[offset] += Fn[t_offset];
3150 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3152 if (negatedFluxNormal[2 * n + 1])
3154 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3159 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3167 static int sav_ncoeffs = 0;
3168 if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3181 m_Ixm = BASE->GetI(coords);
3184 m_Ixp = BASE->GetI(coords);
3186 sav_ncoeffs = e_ncoeffs;
3190 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3192 if (negatedFluxNormal[2 * n])
3194 for (j = 0; j < e_ncoeffs; j++)
3196 outarray[offset + j] -=
3197 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3202 for (j = 0; j < e_ncoeffs; j++)
3204 outarray[offset + j] +=
3205 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3210 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3212 if (negatedFluxNormal[2 * n + 1])
3214 for (j = 0; j < e_ncoeffs; j++)
3216 outarray[offset + j] -=
3217 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3222 for (j = 0; j < e_ncoeffs; j++)
3224 outarray[offset + j] +=
3225 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3239 m_trace->IProductWRTBase(Fn, Fcoeffs);
3246 int e, n, offset, t_offset;
3256 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3258 t_offset =
GetTrace()->GetPhys_Offset(
3259 elmtToTrace[n][e]->GetElmtId());
3260 (*m_exp)[n]->AddEdgeNormBoundaryInt(
3261 e, elmtToTrace[n][e], Fn + t_offset,
3262 e_outarray = outarray + offset);
3271 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3273 t_offset =
GetTrace()->GetPhys_Offset(
3274 elmtToTrace[n][e]->GetElmtId());
3275 (*m_exp)[n]->AddFaceNormBoundaryInt(
3276 e, elmtToTrace[n][e], Fn + t_offset,
3277 e_outarray = outarray + offset);
3311 "tested for 1D expansion");
3315 m_trace->IProductWRTBase(Fwd, Coeffs);
3317 m_trace->IProductWRTBase(Bwd, Coeffs);
3327 const bool PhysSpaceForcing)
3329 boost::ignore_unused(varfactors, dirForcing);
3330 int i, n, cnt, nbndry;
3340 if (PhysSpaceForcing)
3353 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
3354 int NumDirBCs =
m_traceMap->GetNumLocalDirBndCoeffs();
3363 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
3373 for (cnt = n = 0; n < nexp; ++n)
3375 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3377 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3384 Floc =
Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3390 m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3407 locid = bndCondMap[cnt + j];
3408 loclambda[locid] = Sign[locid] * bndcoeffs[j];
3419 locid = bndCondMap[cnt + j];
3420 bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3431 if (GloBndDofs - NumDirBCs > 0)
3436 LinSys->Solve(bndrhs, loclambda,
m_traceMap);
3454 out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3493 const std::string varName,
3509 npoints = locExpList->GetNpoints();
3514 locExpList->GetCoords(x0, x1, x2);
3537 0, (std::static_pointer_cast<
3538 SpatialDomains ::DirichletBoundaryCondition>(
3540 ->m_dirichletCondition)
3541 .Evaluate(x0[0], x1[0], x2[0], time));
3549 0, (std::static_pointer_cast<
3550 SpatialDomains ::NeumannBoundaryCondition>(
3552 ->m_neumannCondition)
3553 .Evaluate(x0[0], x1[0], x2[0], time));
3559 0, (std::static_pointer_cast<
3560 SpatialDomains ::RobinBoundaryCondition>(
3563 .Evaluate(x0[0], x1[0], x2[0], time));
3577 "This type of BC not implemented yet");
3586 std::static_pointer_cast<
3591 valuesExp(npoints, 1.0);
3593 string filebcs = bcPtr->m_filename;
3594 string exprbcs = bcPtr->m_expr;
3600 valuesFile = locExpList->GetPhys();
3606 std::static_pointer_cast<
3609 ->m_dirichletCondition;
3611 condition.
Evaluate(x0, x1, x2, time, valuesExp);
3615 locExpList->UpdatePhys(), 1);
3617 locExpList->FwdTransBndConstrained(
3618 locExpList->GetPhys(), locExpList->UpdateCoeffs());
3624 std::static_pointer_cast<
3636 std::static_pointer_cast<
3639 ->m_neumannCondition;
3640 condition.
Evaluate(x0, x1, x2, time,
3641 locExpList->UpdatePhys());
3644 locExpList->IProductWRTBase(locExpList->GetPhys(),
3645 locExpList->UpdateCoeffs());
3651 std::static_pointer_cast<
3665 std::static_pointer_cast<
3669 condition.
Evaluate(x0, x1, x2, time,
3670 locExpList->UpdatePhys());
3673 locExpList->IProductWRTBase(locExpList->GetPhys(),
3674 locExpList->UpdateCoeffs());
3684 "This type of BC not implemented yet");
3713 id2 =
m_trace->GetPhys_Offset(
3714 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3729 id2 =
m_trace->GetPhys_Offset(
3730 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3741 "Method not set up for this boundary condition.");
3758 map<int, int> VertGID;
3776 VertGID[Vid] = cnt++;
3784 for (i = 0; i < exp->GetNverts(); ++i)
3786 id = exp->GetGeom()->GetVid(i);
3788 if (VertGID.count(
id) > 0)
3790 bid = VertGID.find(
id)->second;
3792 "Edge already set");
3800 "Failed to visit all boundary condtiions");
3805 map<int, int> globalIdMap;
3814 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3840 m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3843 globalIdMap[(*tmp)[0].first->GetGlobalID()];
3851 map<int, int> globalIdMap;
3861 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3885 m_graph->GetElementsFromFace(exp2d->GetGeom2D());
3887 globalIdMap[tmp->at(0).first->GetGlobalID()];
3894 ASSERTL1(
false,
"Not setup for this expansion");
3904 std::shared_ptr<ExpList> &result,
3905 const bool DeclareCoeffPhysArrays)
3908 int offsetOld, offsetNew;
3909 std::vector<unsigned int> eIDs;
3915 for (cnt = n = 0; n < i; ++n)
3923 eIDs.push_back(ElmtID[cnt + n]);
3931 if (DeclareCoeffPhysArrays)
3934 for (n = 0; n < result->GetExpSize(); ++n)
3936 nq =
GetExp(ElmtID[cnt + n])->GetTotPoints();
3938 offsetNew = result->GetPhys_Offset(n);
3940 tmp2 = result->UpdatePhys() + offsetNew, 1);
3942 nq =
GetExp(ElmtID[cnt + n])->GetNcoeffs();
3944 offsetNew = result->GetCoeff_Offset(n);
3946 tmp2 = result->UpdateCoeffs() + offsetNew, 1);
3982 map<int, RobinBCInfoSharedPtr> returnval;
3998 int npoints = locExpList->GetNpoints();
4004 locExpList->GetCoords(x0, x1, x2);
4007 std::static_pointer_cast<
4009 ->m_robinPrimitiveCoeff;
4012 coeffeqn.
Evaluate(x0, x1, x2, 0.0, coeffphys);
4014 for (e = 0; e < locExpList->GetExpSize(); ++e)
4019 Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4021 elmtid = ElmtID[cnt + e];
4023 if (returnval.count(elmtid) != 0)
4025 rInfo->next = returnval.find(elmtid)->second;
4027 returnval[elmtid] = rInfo;
4049 int i, e, ncoeff_edge;
4059 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4064 edge_lambda = loc_lambda;
4072 int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4075 for (e = 0; e < nEdges; ++e)
4077 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4078 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4080 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4081 elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4083 edge_lambda = edge_lambda + ncoeff_edge;
4087 elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4088 cnt += (*m_exp)[i]->GetNcoeffs();
4118 int i, cnt, e, ncoeff_trace;
4125 int nq_elmt, nm_elmt;
4126 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
4131 trace_lambda = loc_lambda;
4135 int num_points[] = {0, 0, 0};
4136 int num_modes[] = {0, 0, 0};
4141 nq_elmt = (*m_exp)[i]->GetTotPoints();
4142 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4146 out_tmp = force + nm_elmt;
4149 for (
int j = 0; j < dim; ++j)
4151 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4152 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4158 int nTraces = (*m_exp)[i]->GetNtraces();
4161 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4163 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4164 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4166 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4169 elmtToTrace[i][e]->SetCoeffsToOrientation(
4170 edgedir, traceCoeffs[e], traceCoeffs[e]);
4176 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4178 trace_lambda = trace_lambda + ncoeff_trace;
4193 num_modes[0], PkeyQ1);
4195 num_modes[1], PkeyQ2);
4197 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4198 (*
m_exp)[i]->GetGeom());
4200 BkeyQ1, BkeyQ2, qGeom);
4208 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4210 num_modes[0], PkeyT1);
4212 num_modes[1], PkeyT2);
4214 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4215 (*
m_exp)[i]->GetGeom());
4217 BkeyT1, BkeyT2, tGeom);
4229 num_modes[0], PkeyH1);
4231 num_modes[1], PkeyH2);
4233 num_modes[2], PkeyH3);
4235 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4236 (*
m_exp)[i]->GetGeom());
4238 BkeyH1, BkeyH2, BkeyH3, hGeom);
4246 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4248 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4250 num_modes[0], PkeyT1);
4252 num_modes[1], PkeyT2);
4254 num_modes[2], PkeyT3);
4256 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4257 (*
m_exp)[i]->GetGeom());
4259 BkeyT1, BkeyT2, BkeyT3, tGeom);
4264 "Wrong shape type, HDG postprocessing is not "
4271 elmtToTrace[i], traceCoeffs, out_tmp);
4272 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4273 ppExp->IProductWRTDerivBase(0, qrhs, force);
4277 elmtToTrace[i], traceCoeffs, out_tmp);
4279 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4280 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4282 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4286 force[0] = (*m_exp)[i]->Integral(qrhs);
4291 ppExp->DetShapeType(), *ppExp);
4297 out = (*lapsys) * in;
4301 ppExp->BwdTrans(out.
GetPtr(), work);
4302 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray +
m_coeff_offset[i]);
4331 m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4333 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...
std::vector< int > m_periodicBwdCopy
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
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
virtual ~DisContField()
Destructor.
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.
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
Array< OneD, int > m_BCtoTraceMap
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
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.
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill the weight with m_bndCondBndWeight.
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
ExpListSharedPtr m_trace
Trace space storage for points between elements.
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
Evaluate all boundary conditions at a given time..
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID)
std::vector< bool > & GetNegatedFluxNormal(void)
DisContField()
Default constructor.
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
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...
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
std::vector< bool > m_leftAdjacentTraces
Array< OneD, int > m_BCtoElmMap
virtual void 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)
Solve the Helmholtz equation.
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.
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
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 global list of m_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.
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
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 with respect to soln of the global spectral/hp element approximat...
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 global list of m_phys correspoinding to element n.
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 BoundaryConditionCollection & GetBoundaryConditions(void) const
const BoundaryRegionCollection & GetBoundaryRegions(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
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.
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::vector< PointGeomSharedPtr > PointGeomVector
std::map< int, CompositeSharedPtr > CompositeMap
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
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< HexGeom > HexGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
std::shared_ptr< Composite > CompositeSharedPtr
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
@ eInvLaplacianWithUnityMean
std::map< ConstFactorType, NekDouble > ConstFactorMap
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1FwdDir1_Dir2FwdDir2
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
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.