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,
84 ExpList3D (pSession, graph3D, variable, ImpType),
85 m_bndCondExpansions(),
90 if (variable.compare(
"DefaultVar") != 0)
117 for(f = 0; f < locExpList->GetExpSize(); ++f)
120 = (*m_exp)[ElmtID[cnt+f]]->
121 as<LocalRegions::Expansion3D>();
123 = locExpList->GetExp(f)->
124 as<LocalRegions::Expansion2D>();
126 exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
127 exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
141 const std::string &variable,
142 const bool SetUpJustDG)
172 for(f = 0; f < locExpList->GetExpSize(); ++f)
175 = (*m_exp)[ElmtID[cnt+f]]->
176 as<LocalRegions::Expansion3D>();
178 = locExpList->GetExp(f)->
179 as<LocalRegions::Expansion2D>();
181 exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
182 exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
217 for(f = 0; f < locExpList->GetExpSize(); ++f)
220 = (*m_exp)[ElmtID[cnt+f]]->
221 as<LocalRegions::Expansion3D>();
223 = locExpList->GetExp(f)->
224 as<LocalRegions::Expansion2D>();
226 exp3d->SetFaceExp(FaceID[cnt+f], exp2d);
227 exp2d->SetAdjacentElementExp(FaceID[cnt+f], exp3d);
233 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
235 std::string ProjectStr =
237 if (ProjectStr ==
"MixedCGDG" ||
238 ProjectStr ==
"Mixed_CG_Discontinuous")
260 m_bndCondExpansions (In.m_bndCondExpansions),
261 m_bndConditions (In.m_bndConditions),
262 m_globalBndMat (In.m_globalBndMat),
263 m_trace (In.m_trace),
264 m_traceMap (In.m_traceMap),
265 m_locTraceToTraceMap (In.m_locTraceToTraceMap),
266 m_periodicFaces (In.m_periodicFaces),
267 m_periodicEdges (In.m_periodicEdges),
268 m_periodicVerts (In.m_periodicVerts)
283 "Routine currently only tested for HybridDGHelmholtz");
286 "The local to global map is not set up for the requested "
295 (*m_globalBndMat)[mkey] = glo_matrix;
299 glo_matrix = matrixIter->second;
326 bool UseGenSegExp =
true;
337 if (
m_session->DefinesCmdLineArgument(
"verbose"))
350 for (
int i = 0; i <
m_exp->size(); ++i)
352 for (
int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
358 exp3d->SetFaceExp (j, exp2d);
367 for (
int i = 0; i <
m_trace->GetExpSize(); ++i)
372 int offset =
m_trace->GetPhys_Offset(i);
373 int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
379 if (traceGeomId != min(pIt->second[0].id, traceGeomId))
381 traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
382 traceEl->GetLeftAdjacentElementFace());
385 else if (
m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
387 traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
388 traceEl->GetLeftAdjacentElementFace());
403 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
410 boost::unordered_map<int,pair<int,int> > perFaceToExpMap;
411 boost::unordered_map<int,pair<int,int> >
::iterator it2;
414 for (
int n = 0; n <
m_exp->size(); ++n)
417 for (
int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
420 exp3d->GetGeom3D()->GetFid(e));
424 perFaceToExpMap[it->first] = make_pair(n, e);
432 for (
int i = 0; i <
m_exp->size(); ++i)
434 for (
int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j, ++cnt)
442 for (
int n = 0; n <
m_exp->size(); ++n)
445 for (
int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
447 int faceGeomId = exp3d->
GetGeom3D()->GetFid(e);
448 int offset =
m_trace->GetPhys_Offset(
449 elmtToTrace[n][e]->GetElmtId());
457 it2 = perFaceToExpMap.find(ent.
id);
459 if (it2 == perFaceToExpMap.end())
461 if (
m_session->GetComm()->GetSize() > 1 &&
468 ASSERTL1(
false,
"Periodic edge not found!");
473 "Periodic face in non-forward space?");
475 int offset2 =
m_trace->GetPhys_Offset(
476 elmtToTrace[it2->second.first][it2->second.second]->
481 int nquad1 = elmtToTrace[n][e]->GetNumPoints(0);
482 int nquad2 = elmtToTrace[n][e]->GetNumPoints(1);
484 vector<int> tmpBwd(nquad1*nquad2);
485 vector<int> tmpFwd(nquad1*nquad2);
492 for (
int i = 0; i < nquad2; ++i)
494 for (
int j = 0; j < nquad1; ++j)
496 tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
497 tmpFwd[i*nquad1+j] = offset + j*nquad2+i;
503 for (
int i = 0; i < nquad2; ++i)
505 for (
int j = 0; j < nquad1; ++j)
507 tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
508 tmpFwd[i*nquad1+j] = offset + i*nquad1+j;
519 for (
int i = 0; i < nquad2; ++i)
521 for (
int j = 0; j < nquad1/2; ++j)
523 swap(tmpFwd[i*nquad1 + j],
524 tmpFwd[i*nquad1 + nquad1-j-1]);
535 for (
int j = 0; j < nquad1; ++j)
537 for (
int i = 0; i < nquad2/2; ++i)
539 swap(tmpFwd[i*nquad1 + j],
540 tmpFwd[(nquad2-i-1)*nquad1 + j]);
545 for (
int i = 0; i < nquad1*nquad2; ++i)
571 bool returnval =
true;
591 int vSame = returnval ? 1 : 0;
612 const std::string &variable)
622 SpatialDomains::BoundaryRegionCollection::const_iterator it;
625 for (it = bregions.begin(); it != bregions.end(); ++it)
629 if (boundaryCondition->GetBoundaryConditionType() !=
642 for (it = bregions.begin(); it != bregions.end(); ++it)
645 bconditions, it->first, variable);
647 if(locBCond->GetBoundaryConditionType()
655 if(locBCond->GetBoundaryConditionType() !=
662 m_bndConditions[cnt++] = locBCond;
676 const std::string &variable)
683 = boost::dynamic_pointer_cast<
685 SpatialDomains::BoundaryRegionCollection::const_iterator it;
708 map<int,int> perComps;
709 map<int,vector<int> > allVerts;
710 map<int,SpatialDomains::PointGeomVector> allCoord;
711 map<int,vector<int> > allEdges;
712 map<int,vector<StdRegions::Orientation> > allOrient;
717 int region1ID, region2ID, i, j, k, cnt;
721 for(i = 0; i < (*m_exp).size(); ++i)
723 for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
725 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
729 for(j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
731 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
739 for (it = bregions.begin(); it != bregions.end(); ++it)
742 bconditions, it->first, variable);
744 if (locBCond->GetBoundaryConditionType()
751 region1ID = it->first;
752 region2ID = boost::static_pointer_cast<
754 locBCond)->m_connectedBoundaryRegion;
758 "Boundary region "+boost::lexical_cast<
string>(
759 region1ID)+
" should only contain 1 composite.");
767 if (vComm->GetSize() == 1)
769 cId1 = it->second->begin()->first;
770 cId2 = bregions.find(region2ID)->second->begin()->first;
774 cId1 = bndRegOrder.find(region1ID)->second[0];
775 cId2 = bndRegOrder.find(region2ID)->second[0];
779 vector<unsigned int> tmpOrder;
785 for (i = 0; i < c->size(); ++i)
788 boost::dynamic_pointer_cast<
790 ASSERTL1(faceGeom,
"Unable to cast to shared ptr");
793 int faceId = (*c)[i]->GetGlobalID();
794 locFaces.insert(faceId);
798 if (vComm->GetSize() == 1)
800 tmpOrder.push_back((*c)[i]->GetGlobalID());
805 vector<int> vertList, edgeList;
807 vector<StdRegions::Orientation> orientVec;
808 for (j = 0; j < faceGeom->GetNumVerts(); ++j)
810 vertList .push_back(faceGeom->GetVid (j));
811 edgeList .push_back(faceGeom->GetEid (j));
812 coordVec .push_back(faceGeom->GetVertex(j));
813 orientVec.push_back(faceGeom->GetEorient(j));
816 allVerts [faceId] = vertList;
817 allEdges [faceId] = edgeList;
818 allCoord [faceId] = coordVec;
819 allOrient[faceId] = orientVec;
824 if (vComm->GetSize() == 1)
826 compOrder[it->second->begin()->first] = tmpOrder;
832 if (perComps.count(cId1) == 0)
834 if (perComps.count(cId2) == 0)
836 perComps[cId1] = cId2;
840 std::stringstream ss;
841 ss <<
"Boundary region " << cId2 <<
" should be "
842 <<
"periodic with " << perComps[cId2] <<
" but "
843 <<
"found " << cId1 <<
" instead!";
844 ASSERTL0(perComps[cId2] == cId1, ss.str());
849 std::stringstream ss;
850 ss <<
"Boundary region " << cId1 <<
" should be "
851 <<
"periodic with " << perComps[cId1] <<
" but "
852 <<
"found " << cId2 <<
" instead!";
853 ASSERTL0(perComps[cId1] == cId1, ss.str());
859 int n = vComm->GetSize();
860 int p = vComm->GetRank();
868 facecounts[
p] = locFaces.size();
874 for (i = 1; i < n; ++i)
876 faceoffset[i] = faceoffset[i-1] + facecounts[i-1];
891 for (i = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
893 faceIds [faceoffset[
p] + i ] = *sIt;
894 faceVerts[faceoffset[
p] + i++] = allVerts[*sIt].size();
917 for (i = 0; i < n; ++i)
919 if (facecounts[i] > 0)
922 facecounts[i], faceVerts + faceoffset[i], 1);
933 for (i = 1; i < n; ++i)
935 vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
949 for (cnt = 0, sIt = locFaces.begin();
950 sIt != locFaces.end(); ++sIt)
952 for (j = 0; j < allVerts[*sIt].size(); ++j)
954 int vertId = allVerts[*sIt][j];
955 vertIds[vertoffset[
p] + cnt ] = vertId;
956 vertX [vertoffset[
p] + cnt ] = (*allCoord[*sIt][j])(0);
957 vertY [vertoffset[
p] + cnt ] = (*allCoord[*sIt][j])(1);
958 vertZ [vertoffset[
p] + cnt ] = (*allCoord[*sIt][j])(2);
959 edgeIds[vertoffset[
p] + cnt ] = allEdges [*sIt][j];
960 edgeOrt[vertoffset[
p] + cnt++] = allOrient[*sIt][j];
975 map<int, vector<int> > vertMap;
976 map<int, vector<int> > edgeMap;
977 map<int, SpatialDomains::PointGeomVector> coordMap;
983 map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
984 map<int, pair<int, int> > eIdMap;
986 for (cnt = i = 0; i < totFaces; ++i)
988 vector<int> edges(faceVerts[i]);
989 vector<int> verts(faceVerts[i]);
995 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
997 edges[j] = edgeIds[cnt];
998 verts[j] = vertIds[cnt];
1001 3, verts[j], vertX[cnt], vertY[cnt], vertZ[cnt]);
1002 vCoMap[vertIds[cnt]] = coord[j];
1005 pair<map<int, pair<int, int> >
::iterator,
bool> testIns =
1006 eIdMap.insert(make_pair(
1008 make_pair(vertIds[tmp+j],
1009 vertIds[tmp+((j+1) % faceVerts[i])])));
1011 if (testIns.second ==
false)
1023 swap(testIns.first->second.first,
1024 testIns.first->second.second);
1028 vertMap [faceIds[i]] = verts;
1029 edgeMap [faceIds[i]] = edges;
1030 coordMap[faceIds[i]] = coord;
1037 map<int,int>::const_iterator oIt;
1051 map<int, map<StdRegions::Orientation, vector<int> > > vmap;
1052 map<int, map<StdRegions::Orientation, vector<int> > > emap;
1054 map<StdRegions::Orientation, vector<int> > quadVertMap;
1064 map<StdRegions::Orientation, vector<int> > quadEdgeMap;
1074 map<StdRegions::Orientation, vector<int> > triVertMap;
1078 map<StdRegions::Orientation, vector<int> > triEdgeMap;
1082 vmap[3] = triVertMap;
1083 vmap[4] = quadVertMap;
1084 emap[3] = triEdgeMap;
1085 emap[4] = quadEdgeMap;
1087 map<int,int> allCompPairs;
1092 for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
1095 const int id1 = cIt->first;
1096 const int id2 = cIt->second;
1097 std::string id1s = boost::lexical_cast<
string>(id1);
1098 std::string id2s = boost::lexical_cast<
string>(id2);
1100 if (compMap.count(id1) > 0)
1102 c[0] = compMap[id1];
1105 if (compMap.count(id2) > 0)
1107 c[1] = compMap[id2];
1111 "Neither composite not found on this process!");
1116 map<int,int> compPairs;
1119 "Unable to find composite "+id1s+
" in order map.");
1121 "Unable to find composite "+id2s+
" in order map.");
1122 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1123 "Periodic composites "+id1s+
" and "+id2s+
1124 " should have the same number of elements.");
1125 ASSERTL0(compOrder[id1].size() > 0,
1126 "Periodic composites "+id1s+
" and "+id2s+
1130 for (i = 0; i < compOrder[id1].size(); ++i)
1132 int eId1 = compOrder[id1][i];
1133 int eId2 = compOrder[id2][i];
1135 ASSERTL0(compPairs.count(eId1) == 0,
1139 if (compPairs.count(eId2) != 0)
1141 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1143 compPairs[eId1] = eId2;
1149 for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1151 int ids [2] = {pIt->first, pIt->second};
1152 bool local[2] = {locFaces.count(pIt->first) > 0,
1153 locFaces.count(pIt->second) > 0};
1155 ASSERTL0(coordMap.count(ids[0]) > 0 &&
1156 coordMap.count(ids[1]) > 0,
1157 "Unable to find face in coordinate map");
1159 allCompPairs[pIt->first ] = pIt->second;
1160 allCompPairs[pIt->second] = pIt->first;
1165 = { coordMap[ids[0]], coordMap[ids[1]] };
1167 ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
1168 "Two periodic faces have different number "
1179 for (i = 0; i < 2; ++i)
1187 int other = (i+1) % 2;
1190 if (tmpVec[0].size() == 3)
1193 tmpVec[i], tmpVec[other]);
1198 tmpVec[i], tmpVec[other]);
1208 int nFaceVerts = vertMap[ids[0]].size();
1211 for (i = 0; i < 2; ++i)
1213 int other = (i+1) % 2;
1216 if (tmpVec[0].size() == 3)
1219 tmpVec[i], tmpVec[other]);
1224 tmpVec[i], tmpVec[other]);
1227 if (nFaceVerts == 3)
1232 "Unsupported face orientation for face "+
1233 boost::lexical_cast<string>(ids[i]));
1237 vector<int> per1 = vertMap[ids[i]];
1238 vector<int> per2 = vertMap[ids[other]];
1242 map<int, pair<int, bool> > tmpMap;
1247 for (j = 0; j < nFaceVerts; ++j)
1249 int v = vmap[nFaceVerts][o][j];
1250 tmpMap[per1[j]] = make_pair(
1251 per2[v], locVerts.count(per2[v]) > 0);
1255 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1259 mIt->second.second);
1265 if (perIt == periodicVerts.end())
1269 periodicVerts[mIt->first].push_back(ent2);
1270 perIt = periodicVerts.find(mIt->first);
1277 for (k = 0; k < perIt->second.size(); ++k)
1279 if (perIt->second[k].id == mIt->second.first)
1285 if (k == perIt->second.size())
1287 perIt->second.push_back(ent2);
1295 for (i = 0; i < 2; ++i)
1297 int other = (i+1) % 2;
1299 if (tmpVec[0].size() == 3)
1302 tmpVec[i], tmpVec[other]);
1307 tmpVec[i], tmpVec[other]);
1310 vector<int> per1 = edgeMap[ids[i]];
1311 vector<int> per2 = edgeMap[ids[other]];
1313 map<int, pair<int, bool> > tmpMap;
1316 for (j = 0; j < nFaceVerts; ++j)
1318 int e = emap[nFaceVerts][o][j];
1319 tmpMap[per1[j]] = make_pair(
1320 per2[e], locEdges.count(per2[e]) > 0);
1323 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1329 mIt->second.second);
1333 if (perIt == periodicEdges.end())
1335 periodicEdges[mIt->first].push_back(ent2);
1336 perIt = periodicEdges.find(mIt->first);
1340 for (k = 0; k < perIt->second.size(); ++k)
1342 if (perIt->second[k].id == mIt->second.first)
1348 if (k == perIt->second.size())
1350 perIt->second.push_back(ent2);
1359 pairSizes[
p] = allCompPairs.size();
1367 for (i = 1; i < n; ++i)
1369 pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1375 cnt = pairOffsets[
p];
1377 for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1379 first [cnt ] = pIt->first;
1380 second[cnt++] = pIt->second;
1386 allCompPairs.clear();
1388 for(cnt = 0; cnt < totPairSizes; ++cnt)
1390 allCompPairs[first[cnt]] = second[cnt];
1397 for (cnt = i = 0; i < totFaces; ++i)
1399 int faceId = faceIds[i];
1401 ASSERTL0(allCompPairs.count(faceId) > 0,
1402 "Unable to find matching periodic face.");
1404 int perFaceId = allCompPairs[faceId];
1406 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1408 int vId = vertIds[cnt];
1412 if (perId == periodicVerts.end())
1420 = { coordMap[faceId], coordMap[perFaceId] };
1422 int nFaceVerts = tmpVec[0].size();
1425 tmpVec[0], tmpVec[1]) :
1427 tmpVec[0], tmpVec[1]);
1431 int perVertexId = vertMap[perFaceId][vmap[nFaceVerts][o][j]];
1436 locVerts.count(perVertexId) > 0);
1438 periodicVerts[vId].push_back(ent);
1441 int eId = edgeIds[cnt];
1443 perId = periodicEdges.find(eId);
1445 if (perId == periodicEdges.end())
1453 = { coordMap[faceId], coordMap[perFaceId] };
1455 int nFaceEdges = tmpVec[0].size();
1458 tmpVec[0], tmpVec[1]) :
1460 tmpVec[0], tmpVec[1]);
1464 int perEdgeId = edgeMap[perFaceId][emap[nFaceEdges][o][j]];
1469 locEdges.count(perEdgeId) > 0);
1471 periodicEdges[eId].push_back(ent);
1479 for (perIt = periodicVerts.begin();
1480 perIt != periodicVerts.end(); ++perIt)
1483 for (i = 0; i < perIt->second.size(); ++i)
1486 perIt2 = periodicVerts.find(perIt->second[i].id);
1487 ASSERTL0(perIt2 != periodicVerts.end(),
1488 "Couldn't find periodic vertex.");
1493 for (j = 0; j < perIt2->second.size(); ++j)
1495 if (perIt2->second[j].id == perIt->first)
1500 for (k = 0; k < perIt->second.size(); ++k)
1502 if (perIt2->second[j].id == perIt->second[k].id)
1508 if (k == perIt->second.size())
1510 perIt->second.push_back(perIt2->second[j]);
1516 for (perIt = periodicEdges.begin();
1517 perIt != periodicEdges.end(); ++perIt)
1519 for (i = 0; i < perIt->second.size(); ++i)
1521 perIt2 = periodicEdges.find(perIt->second[i].id);
1522 ASSERTL0(perIt2 != periodicEdges.end(),
1523 "Couldn't find periodic edge.");
1525 for (j = 0; j < perIt2->second.size(); ++j)
1527 if (perIt2->second[j].id == perIt->first)
1532 for (k = 0; k < perIt->second.size(); ++k)
1534 if (perIt2->second[j].id == perIt->second[k].id)
1540 if (k == perIt->second.size())
1542 perIt->second.push_back(perIt2->second[j]);
1549 for (perIt = periodicEdges.begin();
1550 perIt != periodicEdges.end(); perIt++)
1554 = eIdMap.find(perIt->first);
1556 *vCoMap[eIt->second.first],
1557 *vCoMap[eIt->second.second]
1563 for (i = 0; i < perIt->second.size(); ++i)
1565 eIt = eIdMap.find(perIt->second[i].id);
1568 *vCoMap[eIt->second.first],
1569 *vCoMap[eIt->second.second]
1572 NekDouble cx = 0.5*(w[0](0)-v[0](0)+w[1](0)-v[1](0));
1573 NekDouble cy = 0.5*(w[0](1)-v[0](1)+w[1](1)-v[1](1));
1574 NekDouble cz = 0.5*(w[0](2)-v[0](2)+w[1](2)-v[1](2));
1576 int vMap[2] = {-1,-1};
1577 for (j = 0; j < 2; ++j)
1582 for (k = 0; k < 2; ++k)
1588 if (sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z1-z)*(z1-z))
1598 ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
1599 "Unable to align periodic vertices.");
1600 ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
1601 (vMap[1] == 0 || vMap[1] == 1) &&
1602 (vMap[0] != vMap[1]),
1603 "Unable to align periodic vertices.");
1616 for (perIt = periodicVerts.begin();
1617 perIt != periodicVerts.end(); ++perIt)
1619 if (locVerts.count(perIt->first) > 0)
1625 for (perIt = periodicEdges.begin();
1626 perIt != periodicEdges.end(); ++perIt)
1628 if (locEdges.count(perIt->first) > 0)
1640 as<LocalRegions::Expansion2D>();
1642 int offset =
m_trace->GetPhys_Offset(traceEl->GetElmtId());
1645 if (traceEl->GetLeftAdjacentElementFace () == -1 ||
1646 traceEl->GetRightAdjacentElementFace() == -1)
1656 int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
1662 fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1667 GetTraceToUniversalMapUnique(offset) >= 0;
1671 else if (traceEl->GetLeftAdjacentElementFace () != -1 &&
1672 traceEl->GetRightAdjacentElementFace() != -1)
1676 (traceEl->GetLeftAdjacentElementExp().get()) ==
1681 ASSERTL2(
false,
"Unconnected trace element!");
1723 int n, cnt,
npts, e;
1735 GetNFwdLocTracePts();
1751 id2 =
m_trace->GetPhys_Offset(
1752 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1769 id2 =
m_trace->GetPhys_Offset(
1770 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1785 ASSERTL0(
false,
"Method only set up for Dirichlet, Neumann "
1786 "and Robin conditions.");
1810 "Field is not in physical space.");
1820 Vmath::Zero(outarray.num_elements(), outarray, 1);
1828 m_traceMap->UniversalTraceAssemble(outarray);
1856 m_trace->IProductWRTBase(Fn, Fcoeffs);
1891 m_trace->IProductWRTBase(Fwd,Coeffs);
1893 m_trace->IProductWRTBase(Bwd,Coeffs);
1907 map<int,int> globalIdMap;
1922 globalIdMap[exp3d->GetGeom3D()->GetGlobalID()] = i;
1941 as<LocalRegions::Expansion2D>();
1944 graph3D->GetElementsFromFace(exp2d->GetGeom2D());
1947 m_Element->GetGlobalID()];
1948 m_BCtoFaceMap[cnt] = (*tmp)[0]->m_FaceIndx;
1957 boost::shared_ptr<ExpList> &result,
1958 const bool DeclareCoeffPhysArrays)
1961 int offsetOld, offsetNew;
1962 std::vector<unsigned int> eIDs;
1968 for (cnt = n = 0; n < i; ++n)
1976 eIDs.push_back(ElmtID[cnt+n]);
1982 (*
this, eIDs, DeclareCoeffPhysArrays);
1985 if (DeclareCoeffPhysArrays)
1988 for (n = 0; n < result->GetExpSize(); ++n)
1990 nq =
GetExp(ElmtID[cnt+n])->GetTotPoints();
1992 offsetNew = result->GetPhys_Offset(n);
1994 tmp2 = result->UpdatePhys()+ offsetNew, 1);
1996 nq =
GetExp(ElmtID[cnt+n])->GetNcoeffs();
1998 offsetNew = result->GetCoeff_Offset(n);
2000 tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
2029 const bool PhysSpaceForcing)
2031 int i,j,n,cnt,cnt1,nbndry;
2042 if(PhysSpaceForcing)
2055 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
2056 int NumDirichlet =
m_traceMap->GetNumLocalDirBndCoeffs();
2074 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2083 for(cnt = cnt1 = n = 0; n < nexp; ++n)
2085 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
2087 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
2089 e_l = loc_lambda + cnt1;
2096 Floc =
Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
2114 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2123 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2133 if(GloBndDofs - NumDirichlet > 0)
2150 m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
2153 out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
2170 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2175 m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2176 LocLambda = (*HDGHelm) * LocLambda;
2177 m_traceMap->AssembleBnd(loc_lambda,outarray);
2196 map<int, RobinBCInfoSharedPtr> returnval;
2211 int npoints = locExpList->GetNpoints();
2217 locExpList->GetCoords(x0, x1, x2);
2220 boost::static_pointer_cast<
2225 coeffeqn.
Evaluate(x0, x1, x2, 0.0, coeffphys);
2227 for(e = 0; e < locExpList->GetExpSize(); ++e)
2232 Array_tmp = coeffphys +
2233 locExpList->GetPhys_Offset(e));
2235 elmtid = ElmtID[cnt+e];
2237 if(returnval.count(elmtid) != 0)
2239 rInfo->next = returnval.find(elmtid)->second;
2241 returnval[elmtid] = rInfo;
2273 int i,cnt,f,ncoeff_face;
2278 int nq_elmt, nm_elmt;
2279 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2284 face_lambda = loc_lambda;
2293 nm_elmt = (*m_exp)[i]->GetNcoeffs();
2297 out_tmp = force + nm_elmt;
2300 int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
2301 int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
2302 int num_points2 = (*m_exp)[i]->GetBasis(2)->GetNumPoints();
2303 int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
2304 int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
2305 int num_modes2 = (*m_exp)[i]->GetBasis(2)->GetNumModes();
2310 int nFaces = (*m_exp)[i]->GetNfaces();
2312 for(f = 0; f < nFaces; ++f)
2314 ncoeff_face = elmtToTrace[i][f]->GetNcoeffs();
2316 Vmath::Vcopy(ncoeff_face, face_lambda, 1, faceCoeffs[f], 1);
2317 exp->SetFaceToGeomOrientation(f, faceCoeffs[f]);
2318 face_lambda = face_lambda + ncoeff_face;
2362 ASSERTL0(
false,
"Wrong shape type, HDG postprocessing is not implemented");
2368 (*m_exp)[i]->DGDeriv(
2370 elmtToTrace[i], faceCoeffs, out_tmp);
2371 (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2372 ppExp->IProductWRTDerivBase(0,qrhs,force);
2376 (*m_exp)[i]->DGDeriv(
2378 elmtToTrace[i], faceCoeffs, out_tmp);
2379 (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2380 ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2385 (*m_exp)[i]->DGDeriv(
2387 elmtToTrace[i], faceCoeffs, out_tmp);
2388 (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2389 ppExp->IProductWRTDerivBase(2,qrhs,out_tmp);
2393 (*m_exp)[i]->BwdTrans(
2395 force[0] = (*m_exp)[i]->Integral(qrhs);
2409 ppExp->BwdTrans(out.
GetPtr(), work);
2410 (*m_exp)[i]->FwdTrans(work,
2449 const std::string varName,
2458 for (i = 0; i < nbnd; ++i)
2463 npoints = locExpList->GetNpoints();
2470 locExpList->GetCoords(x0, x1, x2);
2479 string exprbcs = bcPtr->m_expr;
2484 valuesFile = locExpList->GetPhys();
2493 condition.
Evaluate(x0, x1, x2, time, valuesExp);
2496 Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1, locExpList->UpdatePhys(), 1);
2498 locExpList->FwdTrans_BndConstrained(
2499 locExpList->GetPhys(),
2500 locExpList->UpdateCoeffs());
2522 condition.
Evaluate(x0, x1, x2, time,
2523 locExpList->UpdatePhys());
2525 locExpList->IProductWRTBase(locExpList->GetPhys(),
2526 locExpList->UpdateCoeffs());
2548 condition.
Evaluate(x0, x1, x2, time,
2549 locExpList->UpdatePhys());
2553 locExpList->IProductWRTBase(locExpList->GetPhys(),
2554 locExpList->UpdateCoeffs());
2559 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)
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
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...
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
GlobalLinSysMapShPtr m_globalBndMat
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
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< NeumannBoundaryCondition > NeumannBCShPtr
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.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
boost::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
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
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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 .
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_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)
Array< OneD, int > m_BCtoElmMap
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, const bool PhysSpaceForcing)
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Array< OneD, int > m_BCtoFaceMap
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()
boost::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
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