46 #include <boost/assign/std/vector.hpp>
47 #include <boost/tuple/tuple.hpp>
54 namespace MultiRegions
66 DisContField3D::DisContField3D() :
68 m_bndCondExpansions (),
81 const std::string &variable,
82 const bool SetUpJustDG)
83 :
ExpList3D (pSession, graph3D, variable),
84 m_bndCondExpansions(),
89 if (variable.compare(
"DefaultVar") != 0)
116 for(f = 0; f < locExpList->GetExpSize(); ++f)
119 = (*m_exp)[ElmtID[cnt+f]]->
120 as<LocalRegions::Expansion3D>();
122 = locExpList->GetExp(f)->
123 as<LocalRegions::Expansion2D>();
125 exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
126 exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
140 const std::string &variable,
141 const bool SetUpJustDG)
171 for(f = 0; f < locExpList->GetExpSize(); ++f)
174 = (*m_exp)[ElmtID[cnt+f]]->
175 as<LocalRegions::Expansion3D>();
177 = locExpList->GetExp(f)->
178 as<LocalRegions::Expansion2D>();
180 exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
181 exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
216 for(f = 0; f < locExpList->GetExpSize(); ++f)
219 = (*m_exp)[ElmtID[cnt+f]]->
220 as<LocalRegions::Expansion3D>();
222 = locExpList->GetExp(f)->
223 as<LocalRegions::Expansion2D>();
225 exp3d->SetFaceExp(FaceID[cnt+f], exp2d);
226 exp2d->SetAdjacentElementExp(FaceID[cnt+f], exp3d);
232 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
234 std::string ProjectStr =
236 if (ProjectStr ==
"MixedCGDG" ||
237 ProjectStr ==
"Mixed_CG_Discontinuous")
259 m_bndCondExpansions (In.m_bndCondExpansions),
260 m_bndConditions (In.m_bndConditions),
261 m_globalBndMat (In.m_globalBndMat),
262 m_trace (In.m_trace),
263 m_traceMap (In.m_traceMap),
264 m_locTraceToTraceMap (In.m_locTraceToTraceMap),
265 m_periodicFaces (In.m_periodicFaces),
266 m_periodicEdges (In.m_periodicEdges),
267 m_periodicVerts (In.m_periodicVerts)
282 "Routine currently only tested for HybridDGHelmholtz");
285 "The local to global map is not set up for the requested "
294 (*m_globalBndMat)[mkey] = glo_matrix;
298 glo_matrix = matrixIter->second;
325 bool UseGenSegExp =
true;
344 for (
int i = 0; i <
m_exp->size(); ++i)
346 for (
int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
352 exp3d->SetFaceExp (j, exp2d);
361 for (
int i = 0; i <
m_trace->GetExpSize(); ++i)
366 int offset =
m_trace->GetPhys_Offset(i);
367 int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
373 if (traceGeomId != min(pIt->second[0].id, traceGeomId))
375 traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
376 traceEl->GetLeftAdjacentElementFace());
379 else if (
m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
381 traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
382 traceEl->GetLeftAdjacentElementFace());
397 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
404 boost::unordered_map<int,pair<int,int> > perFaceToExpMap;
405 boost::unordered_map<int,pair<int,int> >
::iterator it2;
408 for (
int n = 0; n <
m_exp->size(); ++n)
411 for (
int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
414 exp3d->GetGeom3D()->GetFid(e));
418 perFaceToExpMap[it->first] = make_pair(n, e);
426 for (
int i = 0; i <
m_exp->size(); ++i)
428 for (
int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j, ++cnt)
436 for (
int n = 0; n <
m_exp->size(); ++n)
439 for (
int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
441 int faceGeomId = exp3d->
GetGeom3D()->GetFid(e);
442 int offset =
m_trace->GetPhys_Offset(
443 elmtToTrace[n][e]->GetElmtId());
451 it2 = perFaceToExpMap.find(ent.
id);
453 if (it2 == perFaceToExpMap.end())
455 if (
m_session->GetComm()->GetSize() > 1 &&
462 ASSERTL1(
false,
"Periodic edge not found!");
467 "Periodic face in non-forward space?");
469 int offset2 =
m_trace->GetPhys_Offset(
470 elmtToTrace[it2->second.first][it2->second.second]->
475 int nquad1 = elmtToTrace[n][e]->GetNumPoints(0);
476 int nquad2 = elmtToTrace[n][e]->GetNumPoints(1);
478 vector<int> tmpBwd(nquad1*nquad2);
479 vector<int> tmpFwd(nquad1*nquad2);
486 for (
int i = 0; i < nquad2; ++i)
488 for (
int j = 0; j < nquad1; ++j)
490 tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
491 tmpFwd[i*nquad1+j] = offset + j*nquad2+i;
497 for (
int i = 0; i < nquad2; ++i)
499 for (
int j = 0; j < nquad1; ++j)
501 tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
502 tmpFwd[i*nquad1+j] = offset + i*nquad1+j;
513 for (
int i = 0; i < nquad2; ++i)
515 for (
int j = 0; j < nquad1/2; ++j)
517 swap(tmpFwd[i*nquad1 + j],
518 tmpFwd[i*nquad1 + nquad1-j-1]);
529 for (
int j = 0; j < nquad1; ++j)
531 for (
int i = 0; i < nquad2/2; ++i)
533 swap(tmpFwd[i*nquad1 + j],
534 tmpFwd[(nquad2-i-1)*nquad1 + j]);
539 for (
int i = 0; i < nquad1*nquad2; ++i)
565 bool returnval =
true;
585 int vSame = returnval ? 1 : 0;
606 const std::string &variable)
616 SpatialDomains::BoundaryRegionCollection::const_iterator it;
619 for (it = bregions.begin(); it != bregions.end(); ++it)
623 if (boundaryCondition->GetBoundaryConditionType() !=
636 for (it = bregions.begin(); it != bregions.end(); ++it)
639 bconditions, it->first, variable);
641 if(locBCond->GetBoundaryConditionType()
649 if(locBCond->GetBoundaryConditionType() !=
656 m_bndConditions[cnt++] = locBCond;
670 const std::string &variable)
677 = boost::dynamic_pointer_cast<
679 SpatialDomains::BoundaryRegionCollection::const_iterator it;
702 map<int,int> perComps;
703 map<int,vector<int> > allVerts;
704 map<int,SpatialDomains::PointGeomVector> allCoord;
705 map<int,vector<int> > allEdges;
706 map<int,vector<StdRegions::Orientation> > allOrient;
711 int region1ID, region2ID, i, j, k, cnt;
715 for(i = 0; i < (*m_exp).size(); ++i)
717 for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
719 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
723 for(j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
725 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
733 for (it = bregions.begin(); it != bregions.end(); ++it)
736 bconditions, it->first, variable);
738 if (locBCond->GetBoundaryConditionType()
745 region1ID = it->first;
746 region2ID = boost::static_pointer_cast<
748 locBCond)->m_connectedBoundaryRegion;
752 "Boundary region "+boost::lexical_cast<
string>(
753 region1ID)+
" should only contain 1 composite.");
761 if (vComm->GetSize() == 1)
763 cId1 = it->second->begin()->first;
764 cId2 = bregions.find(region2ID)->second->begin()->first;
768 cId1 = bndRegOrder.find(region1ID)->second[0];
769 cId2 = bndRegOrder.find(region2ID)->second[0];
773 vector<unsigned int> tmpOrder;
779 for (i = 0; i < c->size(); ++i)
782 boost::dynamic_pointer_cast<
784 ASSERTL1(faceGeom,
"Unable to cast to shared ptr");
787 int faceId = (*c)[i]->GetGlobalID();
788 locFaces.insert(faceId);
792 if (vComm->GetSize() == 1)
794 tmpOrder.push_back((*c)[i]->GetGlobalID());
799 vector<int> vertList, edgeList;
801 vector<StdRegions::Orientation> orientVec;
802 for (j = 0; j < faceGeom->GetNumVerts(); ++j)
804 vertList .push_back(faceGeom->GetVid (j));
805 edgeList .push_back(faceGeom->GetEid (j));
806 coordVec .push_back(faceGeom->GetVertex(j));
807 orientVec.push_back(faceGeom->GetEorient(j));
810 allVerts [faceId] = vertList;
811 allEdges [faceId] = edgeList;
812 allCoord [faceId] = coordVec;
813 allOrient[faceId] = orientVec;
818 if (vComm->GetSize() == 1)
820 compOrder[it->second->begin()->first] = tmpOrder;
826 if (perComps.count(cId1) == 0)
828 if (perComps.count(cId2) == 0)
830 perComps[cId1] = cId2;
834 std::stringstream ss;
835 ss <<
"Boundary region " << cId2 <<
" should be "
836 <<
"periodic with " << perComps[cId2] <<
" but "
837 <<
"found " << cId1 <<
" instead!";
838 ASSERTL0(perComps[cId2] == cId1, ss.str());
843 std::stringstream ss;
844 ss <<
"Boundary region " << cId1 <<
" should be "
845 <<
"periodic with " << perComps[cId1] <<
" but "
846 <<
"found " << cId2 <<
" instead!";
847 ASSERTL0(perComps[cId1] == cId1, ss.str());
853 int n = vComm->GetSize();
854 int p = vComm->GetRank();
862 facecounts[p] = locFaces.size();
868 for (i = 1; i < n; ++i)
870 faceoffset[i] = faceoffset[i-1] + facecounts[i-1];
885 for (i = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
887 faceIds [faceoffset[p] + i ] = *sIt;
888 faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
911 for (i = 0; i < n; ++i)
913 if (facecounts[i] > 0)
916 facecounts[i], faceVerts + faceoffset[i], 1);
927 for (i = 1; i < n; ++i)
929 vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
943 for (cnt = 0, sIt = locFaces.begin();
944 sIt != locFaces.end(); ++sIt)
946 for (j = 0; j < allVerts[*sIt].size(); ++j)
948 int vertId = allVerts[*sIt][j];
949 vertIds[vertoffset[p] + cnt ] = vertId;
950 vertX [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(0);
951 vertY [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(1);
952 vertZ [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(2);
953 edgeIds[vertoffset[p] + cnt ] = allEdges [*sIt][j];
954 edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
969 map<int, vector<int> > vertMap;
970 map<int, vector<int> > edgeMap;
971 map<int, SpatialDomains::PointGeomVector> coordMap;
977 map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
978 map<int, pair<int, int> > eIdMap;
980 for (cnt = i = 0; i < totFaces; ++i)
982 vector<int> edges(faceVerts[i]);
983 vector<int> verts(faceVerts[i]);
989 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
991 edges[j] = edgeIds[cnt];
992 verts[j] = vertIds[cnt];
995 3, verts[j], vertX[cnt], vertY[cnt], vertZ[cnt]);
996 vCoMap[vertIds[cnt]] = coord[j];
999 pair<map<int, pair<int, int> >
::iterator,
bool> testIns =
1000 eIdMap.insert(make_pair(
1002 make_pair(vertIds[tmp+j],
1003 vertIds[tmp+((j+1) % faceVerts[i])])));
1005 if (testIns.second ==
false)
1017 swap(testIns.first->second.first,
1018 testIns.first->second.second);
1022 vertMap [faceIds[i]] = verts;
1023 edgeMap [faceIds[i]] = edges;
1024 coordMap[faceIds[i]] = coord;
1031 map<int,int>::const_iterator oIt;
1045 map<int, map<StdRegions::Orientation, vector<int> > > vmap;
1046 map<int, map<StdRegions::Orientation, vector<int> > > emap;
1048 map<StdRegions::Orientation, vector<int> > quadVertMap;
1058 map<StdRegions::Orientation, vector<int> > quadEdgeMap;
1068 map<StdRegions::Orientation, vector<int> > triVertMap;
1072 map<StdRegions::Orientation, vector<int> > triEdgeMap;
1076 vmap[3] = triVertMap;
1077 vmap[4] = quadVertMap;
1078 emap[3] = triEdgeMap;
1079 emap[4] = quadEdgeMap;
1081 map<int,int> allCompPairs;
1086 for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
1089 const int id1 = cIt->first;
1090 const int id2 = cIt->second;
1091 std::string id1s = boost::lexical_cast<
string>(id1);
1092 std::string id2s = boost::lexical_cast<
string>(id2);
1094 if (compMap.count(id1) > 0)
1096 c[0] = compMap[id1];
1099 if (compMap.count(id2) > 0)
1101 c[1] = compMap[id2];
1105 "Neither composite not found on this process!");
1110 map<int,int> compPairs;
1113 "Unable to find composite "+id1s+
" in order map.");
1115 "Unable to find composite "+id2s+
" in order map.");
1116 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1117 "Periodic composites "+id1s+
" and "+id2s+
1118 " should have the same number of elements.");
1119 ASSERTL0(compOrder[id1].size() > 0,
1120 "Periodic composites "+id1s+
" and "+id2s+
1124 for (i = 0; i < compOrder[id1].size(); ++i)
1126 int eId1 = compOrder[id1][i];
1127 int eId2 = compOrder[id2][i];
1129 ASSERTL0(compPairs.count(eId1) == 0,
1133 if (compPairs.count(eId2) != 0)
1135 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1137 compPairs[eId1] = eId2;
1143 for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1145 int ids [2] = {pIt->first, pIt->second};
1146 bool local[2] = {locFaces.count(pIt->first) > 0,
1147 locFaces.count(pIt->second) > 0};
1149 ASSERTL0(coordMap.count(ids[0]) > 0 &&
1150 coordMap.count(ids[1]) > 0,
1151 "Unable to find face in coordinate map");
1153 allCompPairs[pIt->first ] = pIt->second;
1154 allCompPairs[pIt->second] = pIt->first;
1159 = { coordMap[ids[0]], coordMap[ids[1]] };
1161 ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
1162 "Two periodic faces have different number "
1173 for (i = 0; i < 2; ++i)
1181 int other = (i+1) % 2;
1184 if (tmpVec[0].size() == 3)
1187 tmpVec[i], tmpVec[other]);
1192 tmpVec[i], tmpVec[other]);
1202 int nFaceVerts = vertMap[ids[0]].size();
1205 for (i = 0; i < 2; ++i)
1207 int other = (i+1) % 2;
1210 if (tmpVec[0].size() == 3)
1213 tmpVec[i], tmpVec[other]);
1218 tmpVec[i], tmpVec[other]);
1221 if (nFaceVerts == 3)
1226 "Unsupported face orientation for face "+
1227 boost::lexical_cast<string>(ids[i]));
1231 vector<int> per1 = vertMap[ids[i]];
1232 vector<int> per2 = vertMap[ids[other]];
1236 map<int, pair<int, bool> > tmpMap;
1241 for (j = 0; j < nFaceVerts; ++j)
1243 int v = vmap[nFaceVerts][o][j];
1244 tmpMap[per1[j]] = make_pair(
1245 per2[v], locVerts.count(per2[v]) > 0);
1249 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1253 mIt->second.second);
1259 if (perIt == periodicVerts.end())
1263 periodicVerts[mIt->first].push_back(ent2);
1264 perIt = periodicVerts.find(mIt->first);
1271 for (k = 0; k < perIt->second.size(); ++k)
1273 if (perIt->second[k].id == mIt->second.first)
1279 if (k == perIt->second.size())
1281 perIt->second.push_back(ent2);
1289 for (i = 0; i < 2; ++i)
1291 int other = (i+1) % 2;
1293 if (tmpVec[0].size() == 3)
1296 tmpVec[i], tmpVec[other]);
1301 tmpVec[i], tmpVec[other]);
1304 vector<int> per1 = edgeMap[ids[i]];
1305 vector<int> per2 = edgeMap[ids[other]];
1307 map<int, pair<int, bool> > tmpMap;
1310 for (j = 0; j < nFaceVerts; ++j)
1312 int e = emap[nFaceVerts][o][j];
1313 tmpMap[per1[j]] = make_pair(
1314 per2[e], locEdges.count(per2[e]) > 0);
1317 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1323 mIt->second.second);
1327 if (perIt == periodicEdges.end())
1329 periodicEdges[mIt->first].push_back(ent2);
1330 perIt = periodicEdges.find(mIt->first);
1334 for (k = 0; k < perIt->second.size(); ++k)
1336 if (perIt->second[k].id == mIt->second.first)
1342 if (k == perIt->second.size())
1344 perIt->second.push_back(ent2);
1353 pairSizes[p] = allCompPairs.size();
1361 for (i = 1; i < n; ++i)
1363 pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1369 cnt = pairOffsets[p];
1371 for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1373 first [cnt ] = pIt->first;
1374 second[cnt++] = pIt->second;
1380 allCompPairs.clear();
1382 for(cnt = 0; cnt < totPairSizes; ++cnt)
1384 allCompPairs[first[cnt]] = second[cnt];
1391 for (cnt = i = 0; i < totFaces; ++i)
1393 int faceId = faceIds[i];
1395 ASSERTL0(allCompPairs.count(faceId) > 0,
1396 "Unable to find matching periodic face.");
1398 int perFaceId = allCompPairs[faceId];
1400 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1402 int vId = vertIds[cnt];
1406 if (perId == periodicVerts.end())
1414 = { coordMap[faceId], coordMap[perFaceId] };
1416 int nFaceVerts = tmpVec[0].size();
1419 tmpVec[0], tmpVec[1]) :
1421 tmpVec[0], tmpVec[1]);
1425 int perVertexId = vertMap[perFaceId][vmap[nFaceVerts][o][j]];
1430 locVerts.count(perVertexId) > 0);
1432 periodicVerts[vId].push_back(ent);
1435 int eId = edgeIds[cnt];
1437 perId = periodicEdges.find(eId);
1439 if (perId == periodicEdges.end())
1447 = { coordMap[faceId], coordMap[perFaceId] };
1449 int nFaceEdges = tmpVec[0].size();
1452 tmpVec[0], tmpVec[1]) :
1454 tmpVec[0], tmpVec[1]);
1458 int perEdgeId = edgeMap[perFaceId][emap[nFaceEdges][o][j]];
1463 locEdges.count(perEdgeId) > 0);
1465 periodicEdges[eId].push_back(ent);
1473 for (perIt = periodicVerts.begin();
1474 perIt != periodicVerts.end(); ++perIt)
1477 for (i = 0; i < perIt->second.size(); ++i)
1480 perIt2 = periodicVerts.find(perIt->second[i].id);
1481 ASSERTL0(perIt2 != periodicVerts.end(),
1482 "Couldn't find periodic vertex.");
1487 for (j = 0; j < perIt2->second.size(); ++j)
1489 if (perIt2->second[j].id == perIt->first)
1494 for (k = 0; k < perIt->second.size(); ++k)
1496 if (perIt2->second[j].id == perIt->second[k].id)
1502 if (k == perIt->second.size())
1504 perIt->second.push_back(perIt2->second[j]);
1510 for (perIt = periodicEdges.begin();
1511 perIt != periodicEdges.end(); ++perIt)
1513 for (i = 0; i < perIt->second.size(); ++i)
1515 perIt2 = periodicEdges.find(perIt->second[i].id);
1516 ASSERTL0(perIt2 != periodicEdges.end(),
1517 "Couldn't find periodic edge.");
1519 for (j = 0; j < perIt2->second.size(); ++j)
1521 if (perIt2->second[j].id == perIt->first)
1526 for (k = 0; k < perIt->second.size(); ++k)
1528 if (perIt2->second[j].id == perIt->second[k].id)
1534 if (k == perIt->second.size())
1536 perIt->second.push_back(perIt2->second[j]);
1543 for (perIt = periodicEdges.begin();
1544 perIt != periodicEdges.end(); perIt++)
1548 = eIdMap.find(perIt->first);
1550 *vCoMap[eIt->second.first],
1551 *vCoMap[eIt->second.second]
1557 for (i = 0; i < perIt->second.size(); ++i)
1559 eIt = eIdMap.find(perIt->second[i].id);
1562 *vCoMap[eIt->second.first],
1563 *vCoMap[eIt->second.second]
1566 NekDouble cx = 0.5*(w[0](0)-v[0](0)+w[1](0)-v[1](0));
1567 NekDouble cy = 0.5*(w[0](1)-v[0](1)+w[1](1)-v[1](1));
1568 NekDouble cz = 0.5*(w[0](2)-v[0](2)+w[1](2)-v[1](2));
1570 int vMap[2] = {-1,-1};
1571 for (j = 0; j < 2; ++j)
1576 for (k = 0; k < 2; ++k)
1582 if (sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z1-z)*(z1-z))
1592 ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
1593 "Unable to align periodic vertices.");
1594 ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
1595 (vMap[1] == 0 || vMap[1] == 1) &&
1596 (vMap[0] != vMap[1]),
1597 "Unable to align periodic vertices.");
1610 for (perIt = periodicVerts.begin();
1611 perIt != periodicVerts.end(); ++perIt)
1613 if (locVerts.count(perIt->first) > 0)
1619 for (perIt = periodicEdges.begin();
1620 perIt != periodicEdges.end(); ++perIt)
1622 if (locEdges.count(perIt->first) > 0)
1634 as<LocalRegions::Expansion2D>();
1636 int offset =
m_trace->GetPhys_Offset(traceEl->GetElmtId());
1639 if (traceEl->GetLeftAdjacentElementFace () == -1 ||
1640 traceEl->GetRightAdjacentElementFace() == -1)
1650 int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
1656 fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1661 GetTraceToUniversalMapUnique(offset) >= 0;
1665 else if (traceEl->GetLeftAdjacentElementFace () != -1 &&
1666 traceEl->GetRightAdjacentElementFace() != -1)
1670 (traceEl->GetLeftAdjacentElementExp().get()) ==
1675 ASSERTL2(
false,
"Unconnected trace element!");
1717 int n, cnt,
npts, e;
1729 GetNFwdLocTracePts();
1745 id2 =
m_trace->GetPhys_Offset(
1746 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1763 id2 =
m_trace->GetPhys_Offset(
1764 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1779 ASSERTL0(
false,
"Method only set up for Dirichlet, Neumann "
1780 "and Robin conditions.");
1804 "Field is not in physical space.");
1814 Vmath::Zero(outarray.num_elements(), outarray, 1);
1822 m_traceMap->UniversalTraceAssemble(outarray);
1850 m_trace->IProductWRTBase(Fn, Fcoeffs);
1885 m_trace->IProductWRTBase(Fwd,Coeffs);
1887 m_trace->IProductWRTBase(Bwd,Coeffs);
1899 map<int,int> globalIdMap;
1914 globalIdMap[exp3d->GetGeom3D()->GetGlobalID()] = i;
1924 if(ElmtID.num_elements() != nbcs)
1929 if(FaceID.num_elements() != nbcs)
1940 as<LocalRegions::Expansion2D>();
1943 graph3D->GetElementsFromFace(exp2d->GetGeom2D());
1945 ElmtID[cnt] = globalIdMap[(*tmp)[0]->
1946 m_Element->GetGlobalID()];
1947 FaceID[cnt] = (*tmp)[0]->m_FaceIndx;
1953 boost::shared_ptr<ExpList> &result)
1956 int offsetOld, offsetNew;
1958 std::vector<unsigned int> eIDs;
1964 for (cnt = n = 0; n < i; ++n)
1972 eIDs.push_back(ElmtID[cnt+n]);
1980 for (n = 0; n < result->GetExpSize(); ++n)
1982 nq =
GetExp(ElmtID[cnt+n])->GetTotPoints();
1984 offsetNew = result->GetPhys_Offset(n);
1986 tmp2 = result->UpdatePhys()+ offsetNew, 1);
1988 nq =
GetExp(ElmtID[cnt+n])->GetNcoeffs();
1990 offsetNew = result->GetCoeff_Offset(n);
1992 tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
2021 int i,j,n,cnt,cnt1,nbndry;
2038 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
2039 int NumDirichlet =
m_traceMap->GetNumLocalDirBndCoeffs();
2057 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2066 for(cnt = cnt1 = n = 0; n < nexp; ++n)
2072 e_l = loc_lambda + cnt1;
2079 Floc =
Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
2097 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2106 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2116 if(GloBndDofs - NumDirichlet > 0)
2133 m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
2136 out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
2153 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2158 m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2159 LocLambda = (*HDGHelm) * LocLambda;
2160 m_traceMap->AssembleBnd(loc_lambda,outarray);
2179 map<int, RobinBCInfoSharedPtr> returnval;
2194 for(e = 0; e < locExpList->GetExpSize(); ++e)
2197 elmtid = ElmtID[cnt+e];
2199 if(returnval.count(elmtid) != 0)
2201 rInfo->next = returnval.find(elmtid)->second;
2203 returnval[elmtid] = rInfo;
2235 int i,cnt,f,ncoeff_face;
2240 int eid,nq_elmt, nm_elmt;
2241 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2246 face_lambda = loc_lambda;
2255 nq_elmt = (*m_exp)[eid]->GetTotPoints();
2256 nm_elmt = (*m_exp)[eid]->GetNcoeffs();
2260 out_tmp = force + nm_elmt;
2263 int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints();
2264 int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints();
2265 int num_points2 = (*m_exp)[eid]->GetBasis(2)->GetNumPoints();
2266 int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes();
2267 int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes();
2268 int num_modes2 = (*m_exp)[eid]->GetBasis(2)->GetNumModes();
2273 int nFaces = (*m_exp)[eid]->GetNfaces();
2275 for(f = 0; f < nFaces; ++f)
2277 ncoeff_face = elmtToTrace[eid][f]->GetNcoeffs();
2279 Vmath::Vcopy(ncoeff_face, face_lambda, 1, faceCoeffs[f], 1);
2280 exp->SetFaceToGeomOrientation(f, faceCoeffs[f]);
2281 face_lambda = face_lambda + ncoeff_face;
2325 ASSERTL0(
false,
"Wrong shape type, HDG postprocessing is not implemented");
2331 (*m_exp)[eid]->DGDeriv(
2333 elmtToTrace[eid], faceCoeffs, out_tmp);
2334 (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2335 ppExp->IProductWRTDerivBase(0,qrhs,force);
2339 (*m_exp)[eid]->DGDeriv(
2341 elmtToTrace[eid], faceCoeffs, out_tmp);
2342 (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2343 ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2348 (*m_exp)[eid]->DGDeriv(
2350 elmtToTrace[eid], faceCoeffs, out_tmp);
2351 (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2352 ppExp->IProductWRTDerivBase(2,qrhs,out_tmp);
2356 (*m_exp)[eid]->BwdTrans(
2358 force[0] = (*m_exp)[eid]->Integral(qrhs);
2372 ppExp->BwdTrans(out.
GetPtr(), work);
2373 (*m_exp)[eid]->FwdTrans(work,
2412 const std::string varName,
2421 for (i = 0; i < nbnd; ++i)
2426 npoints = locExpList->GetNpoints();
2433 locExpList->GetCoords(x0, x1, x2);
2438 string filebcs = boost::static_pointer_cast<
2442 string exprbcs = boost::static_pointer_cast<
2449 valuesFile = locExpList->GetPhys();
2458 condition.
Evaluate(x0, x1, x2, time, valuesExp);
2461 Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1, locExpList->UpdatePhys(), 1);
2463 locExpList->FwdTrans_BndConstrained(
2464 locExpList->GetPhys(),
2465 locExpList->UpdateCoeffs());
2470 string filebcs = boost::static_pointer_cast<
2486 condition.
Evaluate(x0, x1, x2, time,
2487 locExpList->UpdatePhys());
2489 locExpList->IProductWRTBase(locExpList->GetPhys(),
2490 locExpList->UpdateCoeffs());
2497 string filebcs = boost::static_pointer_cast<
2512 condition.
Evaluate(x0, x1, x2, time,
2513 locExpList->UpdatePhys());
2522 locExpList->IProductWRTBase(locExpList->GetPhys(),
2523 locExpList->UpdateCoeffs());
2528 locExpList->UpdatePhys());
2533 ASSERTL0(
false,
"This type of BC not implemented yet");
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
boost::shared_ptr< MeshGraph3D > MeshGraph3DSharedPtr
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
#define ASSERTL0(condition, msg)
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
virtual ~DisContField3D()
Destructor.
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
boost::shared_ptr< ElementFaceVector > ElementFaceVectorSharedPtr
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result)
static ExpListSharedPtr NullExpListSharedPtr
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::vector< PointGeomSharedPtr > PointGeomVector
void SetAdjacentElementExp(int face, Expansion3DSharedPtr &f)
bool isLocal
Flag specifying if this entity is local to this partition.
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
bool IsLeftAdjacentFace(const int n, const int e)
const BoundaryConditionCollection & GetBoundaryConditions(void) const
void ExtractFileBCs(const std::string &fileName, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
GlobalLinSysMapShPtr m_globalBndMat
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
std::vector< bool > m_leftAdjacentFaces
std::set< int > m_boundaryFaces
A set storing the global IDs of any boundary faces.
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
std::map< int, std::vector< unsigned int > > BndRegionOrdering
std::map< ConstFactorType, NekDouble > ConstFactorMap
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
std::map< int, std::vector< unsigned int > > CompositeOrdering
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.
void ApplyGeomInfo()
Apply geometry information to each expansion.
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
int GetExpSize(void)
This function returns the number of elements in the expansion.
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)
This function evaluates the boundary conditions at a certain time-level.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
boost::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
Principle Orthogonal Functions .
void FindPeriodicFaces(const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
Determine the periodic faces, edges and vertices for the given graph.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2)
Get the orientation of face1.
The base class for all shapes.
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
std::vector< int > m_periodicBwdCopy
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
DisContField3D()
Default constructor.
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Map of local trace (the points at the face of the element) to the trace space discretisation.
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
NekDouble Evaluate() const
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Principle Orthogonal Functions .
Principle Orthogonal Functions .
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph3D, const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Calculates the result of the multiplication of a global matrix of type specified by mkey with a vecto...
LibUtilities::SessionReaderSharedPtr m_session
Session.
Defines a specification for a set of points.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
boost::shared_ptr< Expansion > ExpansionSharedPtr
bool SameTypeOfBoundaryConditions(const DisContField3D &In)
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
boost::shared_ptr< GeometryVector > Composite
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Describe a linear system.
std::map< int, Composite > CompositeMap
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::CommSharedPtr m_comm
Communicator.
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const Array< OneD, const NekDouble > &dirForcing)
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &FaceID)
Set up a list of elemeent IDs and edge IDs that link to the boundary conditions.
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Gauss Radau pinned at x=-1, .
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
AssemblyMapDGSharedPtr m_traceMap
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Abstraction of a three-dimensional multi-elemental expansion which is merely a collection of local ex...
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2)
const BoundaryRegionCollection & GetBoundaryRegions(void) const
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Describes the specification for a Basis.
Array< OneD, DataType > & GetPtr()
1D Gauss-Lobatto-Legendre quadrature points
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
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.
PeriodicMap m_periodicEdges
A map which identifies groups of periodic edges.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr