51 namespace MultiRegions
63 DisContField2D::DisContField2D(
void)
65 m_bndCondExpansions(),
73 const bool DeclareCoeffPhysArrays)
75 m_bndCondExpansions (In.m_bndCondExpansions),
76 m_bndConditions (In.m_bndConditions),
77 m_globalBndMat (In.m_globalBndMat),
78 m_traceMap (In.m_traceMap),
79 m_locTraceToTraceMap (In.m_locTraceToTraceMap),
80 m_boundaryEdges (In.m_boundaryEdges),
81 m_periodicVerts (In.m_periodicVerts),
82 m_periodicEdges (In.m_periodicEdges),
83 m_periodicFwdCopy (In.m_periodicFwdCopy),
84 m_periodicBwdCopy (In.m_periodicBwdCopy),
85 m_leftAdjacentEdges (In.m_leftAdjacentEdges)
90 *boost::dynamic_pointer_cast<ExpList1D>(In.
m_trace),
91 DeclareCoeffPhysArrays);
102 const std::string &variable,
103 const bool SetUpJustDG,
104 const bool DeclareCoeffPhysArrays,
106 :
ExpList2D(pSession, graph2D, DeclareCoeffPhysArrays, variable,
108 m_bndCondExpansions(),
117 if (variable.compare(
"DefaultVar") != 0)
121 DeclareCoeffPhysArrays);
123 if (DeclareCoeffPhysArrays)
149 for(e = 0; e < locExpList->GetExpSize(); ++e)
152 (*m_exp)[ElmtID[cnt+e]]->
153 as<LocalRegions::Expansion2D>();
155 locExpList->GetExp(e)->
156 as<LocalRegions::Expansion1D>();
158 locExpList->GetExp(e)->
159 as<LocalRegions::Expansion> ();
161 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
162 exp1d->SetAdjacentElementExp(EdgeID[cnt+e], exp2d);
167 if(
m_session->DefinesSolverInfo(
"PROJECTION"))
169 std::string ProjectStr =
171 if((ProjectStr ==
"MixedCGDG") ||
172 (ProjectStr ==
"Mixed_CG_Discontinuous"))
196 const std::string &variable,
197 const bool SetUpJustDG,
198 const bool DeclareCoeffPhysArrays):
204 if(variable.compare(
"DefaultVar") != 0)
209 if (DeclareCoeffPhysArrays)
237 for(e = 0; e < locExpList->GetExpSize(); ++e)
240 = (*m_exp)[ElmtID[cnt+e]]->
241 as<LocalRegions::Expansion2D>();
243 = locExpList->GetExp(e)->
244 as<LocalRegions::Expansion1D>();
246 = locExpList->GetExp(e)->
247 as<LocalRegions::Expansion> ();
249 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
250 exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
258 if (
m_session->DefinesSolverInfo(
"PROJECTION"))
260 std::string ProjectStr =
263 if ((ProjectStr ==
"MixedCGDG") ||
264 (ProjectStr ==
"Mixed_CG_Discontinuous"))
319 for (e = 0; e < locExpList->GetExpSize(); ++e)
322 = (*m_exp)[ElmtID[cnt+e]]->
323 as<LocalRegions::Expansion2D>();
325 = locExpList->GetExp(e)->
326 as<LocalRegions::Expansion1D>();
328 = locExpList->GetExp(e)->
329 as<LocalRegions::Expansion> ();
331 exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
332 exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
355 "Routine currently only tested for HybridDGHelmholtz");
358 "The local to global map is not set up for the requested "
367 (*m_globalBndMat)[mkey] = glo_matrix;
371 glo_matrix = matrixIter->second;
409 if (
m_session->DefinesCmdLineArgument(
"verbose"))
411 m_traceMap->PrintStats(std::cout, variable);
415 &elmtToTrace = m_traceMap->GetElmtToTrace();
422 for (
int i = 0; i <
m_exp->size(); ++i)
424 for (
int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
431 exp2d->SetEdgeExp (j, exp1d);
432 exp1d->SetAdjacentElementExp(j, exp2d);
440 for (
int i = 0; i <
m_trace->GetExpSize(); ++i)
445 int offset =
m_trace->GetPhys_Offset(i);
446 int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
451 if (traceGeomId != min(pIt->second[0].id, traceGeomId))
453 traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
454 traceEl->GetLeftAdjacentElementEdge());
457 else if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
459 traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
460 traceEl->GetLeftAdjacentElementEdge());
475 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
482 boost::unordered_map<int,pair<int,int> > perEdgeToExpMap;
483 boost::unordered_map<int,pair<int,int> >
::iterator it2;
484 for (cnt = n = 0; n <
m_exp->size(); ++n)
486 for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
489 (*
m_exp)[n]->GetGeom()->GetEid(e));
493 perEdgeToExpMap[it->first] = make_pair(n, e);
501 for (
int i = 0; i <
m_exp->size(); ++i)
503 for (
int j = 0; j < (*m_exp)[i]->GetNedges(); ++j, ++cnt)
511 for (
int n = 0; n <
m_exp->size(); ++n)
513 for (
int e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
515 int edgeGeomId = (*m_exp)[n]->GetGeom()->GetEid(e);
516 int offset =
m_trace->GetPhys_Offset(
517 elmtToTrace[n][e]->GetElmtId());
525 it2 = perEdgeToExpMap.find(ent.
id);
527 if (it2 == perEdgeToExpMap.end())
530 GetRowComm()->GetSize() > 1 && !ent.
isLocal)
536 ASSERTL1(
false,
"Periodic edge not found!");
541 "Periodic edge in non-forward space?");
543 int offset2 =
m_trace->GetPhys_Offset(
544 elmtToTrace[it2->second.first][it2->second.second]->
549 int nquad = elmtToTrace[n][e]->GetNumPoints(0);
551 vector<int> tmpBwd(nquad);
552 vector<int> tmpFwd(nquad);
556 for (
int i = 0; i < nquad; ++i)
558 tmpBwd[i] = offset2 + i;
559 tmpFwd[i] = offset + i;
564 for (
int i = 0; i < nquad; ++i)
566 tmpBwd[i] = offset2 + i;
567 tmpFwd[i] = offset + nquad - i - 1;
571 for (
int i = 0; i < nquad; ++i)
595 bool returnval =
true;
615 int vSame = (returnval?1:0);
639 const std::string &variable,
640 const bool DeclareCoeffPhysArrays)
649 SpatialDomains::BoundaryRegionCollection::const_iterator it;
652 for (it = bregions.begin(); it != bregions.end(); ++it)
670 for (it = bregions.begin(); it != bregions.end(); ++it)
678 DeclareCoeffPhysArrays, variable);
681 m_bndConditions[cnt] = bc;
684 std::string type = m_bndConditions[cnt]->GetUserDefined();
689 if((bc->GetBoundaryConditionType() !=
691 boost::iequals(type,
"I") ||
692 boost::iequals(type,
"CalcBC"))
715 const std::string &variable)
722 = boost::dynamic_pointer_cast<
724 SpatialDomains::BoundaryRegionCollection::const_iterator it;
738 map<int,int> perComps;
739 map<int,vector<int> > allVerts;
741 map<int,StdRegions::Orientation> allEdges;
743 int region1ID, region2ID, i, j, k, cnt;
747 for(i = 0; i < (*m_exp).size(); ++i)
749 for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
751 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
757 for (it = bregions.begin(); it != bregions.end(); ++it)
760 bconditions, it->first, variable);
762 if (locBCond->GetBoundaryConditionType()
769 region1ID = it->first;
770 region2ID = boost::static_pointer_cast<
772 locBCond)->m_connectedBoundaryRegion;
777 if (vComm->GetSize() == 1)
779 cId1 = it->second->begin()->first;
780 cId2 = bregions.find(region2ID)->second->begin()->first;
784 cId1 = bndRegOrder.find(region1ID)->second[0];
785 cId2 = bndRegOrder.find(region2ID)->second[0];
789 "Boundary region "+boost::lexical_cast<
string>(
790 region1ID)+
" should only contain 1 composite.");
795 vector<unsigned int> tmpOrder;
797 for (i = 0; i < c->size(); ++i)
800 boost::dynamic_pointer_cast<
802 ASSERTL0(segGeom,
"Unable to cast to shared ptr");
805 graph2D->GetElementsFromEdge(segGeom);
807 "The periodic boundaries belong to "
808 "more than one element of the mesh");
812 (*elmt)[0]->m_Element);
814 allEdges[(*c)[i]->GetGlobalID()] =
815 geom->GetEorient((*elmt)[0]->m_EdgeIndx);
819 if (vComm->GetSize() == 1)
821 tmpOrder.push_back((*c)[i]->GetGlobalID());
824 vector<int> vertList(2);
825 vertList[0] = segGeom->GetVid(0);
826 vertList[1] = segGeom->GetVid(1);
827 allVerts[(*c)[i]->GetGlobalID()] = vertList;
830 if (vComm->GetSize() == 1)
832 compOrder[it->second->begin()->first] = tmpOrder;
837 if (perComps.count(cId1) == 0)
839 if (perComps.count(cId2) == 0)
841 perComps[cId1] = cId2;
845 std::stringstream ss;
846 ss <<
"Boundary region " << cId2 <<
" should be "
847 <<
"periodic with " << perComps[cId2] <<
" but "
848 <<
"found " << cId1 <<
" instead!";
849 ASSERTL0(perComps[cId2] == cId1, ss.str());
854 std::stringstream ss;
855 ss <<
"Boundary region " << cId1 <<
" should be "
856 <<
"periodic with " << perComps[cId1] <<
" but "
857 <<
"found " << cId2 <<
" instead!";
858 ASSERTL0(perComps[cId1] == cId1, ss.str());
863 int n = vComm->GetSize();
864 int p = vComm->GetRank();
870 edgecounts[
p] = allEdges.size();
874 for (i = 1; i < n; ++i)
876 edgeoffset[i] = edgeoffset[i-1] + edgecounts[i-1];
886 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
888 edgeIds [edgeoffset[
p] + i ] = sIt->first;
889 edgeOrient[edgeoffset[
p] + i ] = sIt->second;
890 edgeVerts [edgeoffset[
p] + i++] = allVerts[sIt->first].size();
912 for (i = 0; i < n; ++i)
914 if (edgecounts[i] > 0)
917 edgecounts[i], edgeVerts + edgeoffset[i], 1);
926 for (i = 1; i < n; ++i)
928 vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
932 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
934 for (j = 0; j < allVerts[sIt->first].size(); ++j)
936 vertIds[vertoffset[
p] + i++] = allVerts[sIt->first][j];
943 map<int, StdRegions::Orientation> orientMap;
944 map<int, vector<int> > vertMap;
946 for (cnt = i = 0; i < totEdges; ++i)
948 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
949 "Already found edge in orientation map!");
952 vector<int> verts(edgeVerts[i]);
954 for (j = 0; j < edgeVerts[i]; ++j)
956 verts[j] = vertIds[cnt++];
958 vertMap[edgeIds[i]] = verts;
965 map<int,int>::const_iterator oIt;
967 map<int,int> allCompPairs;
976 for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
979 const int id1 = cIt->first;
980 const int id2 = cIt->second;
981 std::string id1s = boost::lexical_cast<
string>(id1);
982 std::string id2s = boost::lexical_cast<
string>(id2);
984 if (compMap.count(id1) > 0)
989 if (compMap.count(id2) > 0)
995 "Both composites not found on this process!");
1000 map<int,int> compPairs;
1003 "Unable to find composite "+id1s+
" in order map.");
1005 "Unable to find composite "+id2s+
" in order map.");
1006 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1007 "Periodic composites "+id1s+
" and "+id2s+
1008 " should have the same number of elements.");
1009 ASSERTL0(compOrder[id1].size() > 0,
1010 "Periodic composites "+id1s+
" and "+id2s+
1014 for (i = 0; i < compOrder[id1].size(); ++i)
1016 int eId1 = compOrder[id1][i];
1017 int eId2 = compOrder[id2][i];
1019 ASSERTL0(compPairs.count(eId1) == 0,
1022 if (compPairs.count(eId2) != 0)
1024 ASSERTL0(compPairs[eId2] == eId1,
"Pairing incorrect");
1026 compPairs[eId1] = eId2;
1033 for (i = 0; i < 2; ++i)
1040 if (c[i]->size() > 0)
1042 for (j = 0; j < c[i]->size(); ++j)
1044 locEdges.insert(c[i]->at(j)->GetGlobalID());
1050 for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1052 int ids [2] = {pIt->first, pIt->second};
1053 bool local[2] = {locEdges.count(pIt->first) > 0,
1054 locEdges.count(pIt->second) > 0};
1056 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1057 orientMap.count(ids[1]) > 0,
1058 "Unable to find edge in orientation map.");
1060 allCompPairs[pIt->first ] = pIt->second;
1061 allCompPairs[pIt->second] = pIt->first;
1063 for (i = 0; i < 2; ++i)
1070 int other = (i+1) % 2;
1073 orientMap[ids[i]] == orientMap[ids[other]] ?
1082 for (i = 0; i < 2; ++i)
1084 int other = (i+1) % 2;
1087 orientMap[ids[i]] == orientMap[ids[other]] ?
1092 vector<int> perVerts1 = vertMap[ids[i]];
1093 vector<int> perVerts2 = vertMap[ids[other]];
1095 map<int, pair<int, bool> > tmpMap;
1099 tmpMap[perVerts1[0]] = make_pair(
1100 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1101 tmpMap[perVerts1[1]] = make_pair(
1102 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1106 tmpMap[perVerts1[0]] = make_pair(
1107 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1108 tmpMap[perVerts1[1]] = make_pair(
1109 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1112 for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1117 mIt->second.second);
1121 if (perIt == periodicVerts.end())
1123 periodicVerts[mIt->first].push_back(ent2);
1124 perIt = periodicVerts.find(mIt->first);
1129 for (j = 0; j < perIt->second.size(); ++j)
1131 if (perIt->second[j].id == mIt->second.first)
1140 perIt->second.push_back(ent2);
1149 pairSizes[
p] = allCompPairs.size();
1157 for (i = 1; i < n; ++i)
1159 pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1165 cnt = pairOffsets[
p];
1167 for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1169 first [cnt ] = pIt->first;
1170 second[cnt++] = pIt->second;
1176 allCompPairs.clear();
1178 for(cnt = 0; cnt < totPairSizes; ++cnt)
1180 allCompPairs[first[cnt]] = second[cnt];
1186 for (cnt = i = 0; i < totEdges; ++i)
1188 int edgeId = edgeIds[i];
1190 ASSERTL0(allCompPairs.count(edgeId) > 0,
1191 "Unable to find matching periodic edge.");
1193 int perEdgeId = allCompPairs[edgeId];
1195 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1197 int vId = vertIds[cnt];
1201 if (perId == periodicVerts.end())
1208 orientMap[edgeId] == orientMap[perEdgeId] ?
1209 vertMap[perEdgeId][(j+1)%2] : vertMap[perEdgeId][j];
1213 locVerts.count(perVertexId) > 0);
1215 periodicVerts[vId].push_back(ent);
1223 for (perIt = periodicVerts.begin();
1224 perIt != periodicVerts.end(); ++perIt)
1227 for (i = 0; i < perIt->second.size(); ++i)
1229 perIt2 = periodicVerts.find(perIt->second[i].id);
1230 ASSERTL0(perIt2 != periodicVerts.end(),
1231 "Couldn't find periodic vertex.");
1233 for (j = 0; j < perIt2->second.size(); ++j)
1235 if (perIt2->second[j].id == perIt->first)
1241 for (k = 0; k < perIt->second.size(); ++k)
1243 if (perIt2->second[j].id == perIt->second[k].id)
1252 perIt->second.push_back(perIt2->second[j]);
1260 for (perIt = periodicVerts.begin();
1261 perIt != periodicVerts.end(); ++perIt)
1263 if (locVerts.count(perIt->first) > 0)
1275 as<LocalRegions::Expansion1D>();
1279 if (traceEl->GetLeftAdjacentElementEdge () == -1 ||
1280 traceEl->GetRightAdjacentElementEdge() == -1)
1290 int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
1296 fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1300 int offset =
m_trace->GetPhys_Offset(traceEl->GetElmtId());
1303 GetTraceToUniversalMapUnique(offset) >= 0;
1307 else if (traceEl->GetLeftAdjacentElementEdge () != -1 &&
1308 traceEl->GetRightAdjacentElementEdge() != -1)
1312 (traceEl->GetLeftAdjacentElementExp().get()) ==
1317 ASSERTL2(
false,
"Unconnected trace element!");
1364 int cnt, n, e,
npts, phys_offset;
1383 GetNFwdLocTracePts();
1392 boost::unordered_map<int,pair<int,int> >
::iterator it3;
1398 for(cnt = n = 0; n < nexp; ++n)
1403 for(e = 0; e < exp2d->GetNedges(); ++e, ++cnt)
1405 int offset =
m_trace->GetPhys_Offset(
1406 elmtToTrace[n][e]->GetElmtId());
1410 exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1411 field + phys_offset,
1412 e_tmp = Fwd + offset);
1416 exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1417 field + phys_offset,
1418 e_tmp = Bwd + offset);
1438 GetBndCondTraceToGlobalTraceMap(cnt+e));
1456 "Method not set up for non-zero Neumann "
1457 "boundary condition");
1458 id2 =
m_trace->GetPhys_Offset(
1459 m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1469 "Method not set up for this boundary condition.");
1489 "Field must be in physical state to extract trace space.");
1511 int n, e, offset, phys_offset;
1517 "input array is of insufficient length");
1520 for(n = 0; n < nexp; ++n)
1524 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1526 offset =
m_trace->GetPhys_Offset(
1527 elmtToTrace[n][e]->GetElmtId());
1528 (*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
1529 inarray + phys_offset,
1530 e_tmp = outarray + offset);
1540 int e, n, offset, t_offset;
1548 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1550 t_offset =
GetTrace()->GetPhys_Offset(
1551 elmtToTrace[n][e]->GetElmtId());
1553 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1554 e,elmtToTrace[n][e],
1557 e_outarray = outarray+offset);
1588 m_trace->IProductWRTBase(Fn, Fcoeffs);
1595 int e, n, offset, t_offset;
1603 for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1605 t_offset =
GetTrace()->GetPhys_Offset(
1606 elmtToTrace[n][e]->GetElmtId());
1607 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1608 e, elmtToTrace[n][e],
1610 e_outarray = outarray+offset);
1644 int e,n,offset, t_offset;
1652 for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1654 t_offset =
GetTrace()->GetPhys_Offset(
1655 elmtToTrace[n][e]->GetElmtId());
1660 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1661 e, elmtToTrace[n][e], Fwd+t_offset,
1662 e_outarray = outarray+offset);
1666 (*m_exp)[n]->AddEdgeNormBoundaryInt(
1667 e, elmtToTrace[n][e], Bwd+t_offset,
1668 e_outarray = outarray+offset);
1685 map<int, int> globalIdMap;
1698 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1718 as<LocalRegions::Expansion1D>();
1721 graph2D->GetElementsFromEdge(exp1d->GetGeom1D());
1724 m_Element->GetGlobalID()];
1725 m_BCtoEdgMap[cnt] = (*tmp)[0]->m_EdgeIndx;
1734 boost::shared_ptr<ExpList> &result,
1735 const bool DeclareCoeffPhysArrays)
1738 int offsetOld, offsetNew;
1739 std::vector<unsigned int> eIDs;
1745 for (cnt = n = 0; n < i; ++n)
1753 eIDs.push_back(ElmtID[cnt+n]);
1759 (*
this, eIDs, DeclareCoeffPhysArrays);
1762 if( DeclareCoeffPhysArrays)
1765 for (n = 0; n < result->GetExpSize(); ++n)
1767 nq =
GetExp(ElmtID[cnt+n])->GetTotPoints();
1769 offsetNew = result->GetPhys_Offset(n);
1771 tmp2 = result->UpdatePhys()+ offsetNew, 1);
1773 nq =
GetExp(ElmtID[cnt+n])->GetNcoeffs();
1775 offsetNew = result->GetCoeff_Offset(n);
1777 tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
1810 int i,e,ncoeff_edge;
1820 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
1824 edge_lambda = loc_lambda;
1832 int nEdges = (*m_exp)[i]->GetNedges();
1835 for(e = 0; e < nEdges; ++e)
1837 edgedir = (*m_exp)[i]->GetEorient(e);
1838 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
1840 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
1841 elmtToTrace[i][e]->SetCoeffsToOrientation(
1842 edgedir, edgeCoeffs[e], edgeCoeffs[e]);
1843 edge_lambda = edge_lambda + ncoeff_edge;
1846 (*m_exp)[i]->DGDeriv(dir,
1850 out_tmp = out_d+cnt);
1851 cnt += (*m_exp)[i]->GetNcoeffs();
1866 const bool PhysSpaceForcing)
1869 int i,j,n,cnt,cnt1,nbndry;
1880 if(PhysSpaceForcing)
1893 int GloBndDofs =
m_traceMap->GetNumGlobalBndCoeffs();
1894 int NumDirichlet =
m_traceMap->GetNumLocalDirBndCoeffs();
1915 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
1924 for(cnt = cnt1 = n = 0; n < nexp; ++n)
1926 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1928 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1930 e_l = loc_lambda + cnt1;
1937 Floc =
Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1956 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1965 id =
m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1975 if(GloBndDofs - NumDirichlet > 0)
1994 m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1997 out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
2015 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2020 m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2021 LocLambda = (*HDGHelm) * LocLambda;
2022 m_traceMap->AssembleBnd(loc_lambda,outarray);
2040 map<int, RobinBCInfoSharedPtr> returnval;
2056 int npoints = locExpList->GetNpoints();
2062 locExpList->GetCoords(x0, x1, x2);
2065 boost::static_pointer_cast<
2070 coeffeqn.
Evaluate(x0, x1, x2, 0.0, coeffphys);
2072 for(e = 0; e < locExpList->GetExpSize(); ++e)
2078 Array_tmp = coeffphys +
2079 locExpList->GetPhys_Offset(e));
2081 elmtid = ElmtID[cnt+e];
2083 if(returnval.count(elmtid) != 0)
2085 rInfo->next = returnval.find(elmtid)->second;
2087 returnval[elmtid] = rInfo;
2119 int i,cnt,e,ncoeff_edge;
2126 int nq_elmt, nm_elmt;
2127 int LocBndCoeffs =
m_traceMap->GetNumLocalBndCoeffs();
2132 edge_lambda = loc_lambda;
2137 nq_elmt = (*m_exp)[i]->GetTotPoints();
2138 nm_elmt = (*m_exp)[i]->GetNcoeffs();
2142 out_tmp = force + nm_elmt;
2145 int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
2146 int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
2147 int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
2148 int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
2153 int nEdges = (*m_exp)[i]->GetNedges();
2156 for(e = 0; e < (*m_exp)[i]->GetNedges(); ++e)
2158 edgedir = (*m_exp)[i]->GetEorient(e);
2159 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
2161 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
2162 elmtToTrace[i][e]->SetCoeffsToOrientation(
2163 edgedir, edgeCoeffs[e], edgeCoeffs[e]);
2164 edge_lambda = edge_lambda + ncoeff_edge;
2192 ASSERTL0(
false,
"Wrong shape type, HDG postprocessing is not implemented");
2198 (*m_exp)[i]->DGDeriv(
2200 elmtToTrace[i], edgeCoeffs, out_tmp);
2201 (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2203 ppExp->IProductWRTDerivBase(0,qrhs,force);
2207 (*m_exp)[i]->DGDeriv(
2209 elmtToTrace[i], edgeCoeffs, out_tmp);
2211 (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2213 ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2218 (*m_exp)[i]->BwdTrans(
2220 force[0] = (*m_exp)[i]->Integral(qrhs);
2234 ppExp->BwdTrans(out.
GetPtr(), work);
2235 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray +
m_coeff_offset[i]);
2253 const std::string varName,
2263 for (i = 0; i < nbnd; ++i)
2269 npoints = locExpList->GetNpoints();
2277 locExpList->GetCoords(x0, x1, x2);
2281 locExpList->GetCoords(x0, x1, x2);
2289 = boost::static_pointer_cast<
2301 boost::static_pointer_cast<
2304 m_dirichletCondition;
2306 condition.
Evaluate(x0, x1, x2, time,
2307 locExpList->UpdatePhys());
2310 locExpList->FwdTrans_BndConstrained(
2311 locExpList->GetPhys(),
2312 locExpList->UpdateCoeffs());
2328 boost::static_pointer_cast<
2332 condition.
Evaluate(x0, x1, x2, time,
2333 locExpList->UpdatePhys());
2336 locExpList->IProductWRTBase(
2337 locExpList->GetPhys(),
2338 locExpList->UpdateCoeffs());
2355 boost::static_pointer_cast<
2359 condition.
Evaluate(x0, x1, x2, time,
2360 locExpList->UpdatePhys());
2363 locExpList->IProductWRTBase(
2364 locExpList->GetPhys(),
2365 locExpList->UpdateCoeffs());
2369 ASSERTL0(
false,
"This type of BC not implemented yet");
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
virtual void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This method extracts the trace (edges in 2D) from the field inarray and puts the values in outarray...
#define ASSERTL0(condition, msg)
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
void FindPeriodicEdges(const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
Determine the periodic edges and vertices for the given graph.
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
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.
Array< OneD, int > m_BCtoElmMap
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...
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
const BoundaryConditionCollection & GetBoundaryConditions(void) const
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
Lagrange Polynomials using the Gauss points .
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
std::map< int, std::vector< unsigned int > > BndRegionOrdering
std::map< ConstFactorType, NekDouble > ConstFactorMap
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
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.
std::map< int, std::vector< unsigned int > > CompositeOrdering
boost::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
boost::shared_ptr< MeshGraph2D > MeshGraph2DSharedPtr
int GetExpSize(void)
This function returns the number of elements in the expansion.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
Search through the edge expansions and identify which ones have Robin/Mixed type boundary conditions...
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, .
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Principle Orthogonal Functions .
boost::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
bool SameTypeOfBoundaryConditions(const DisContField2D &In)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Base class for all multi-elemental spectral/hp expansions.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
DisContField2D()
Default constructor.
The base class for all shapes.
boost::shared_ptr< ExpList > & GetTrace()
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
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.
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)
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
virtual ~DisContField2D()
Default destructor.
std::vector< bool > m_leftAdjacentEdges
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
boost::shared_ptr< SegGeom > SegGeomSharedPtr
std::set< int > m_boundaryEdges
A set storing the global IDs of any boundary edges.
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 .
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)
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Principle Orthogonal Functions .
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
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.
boost::shared_ptr< Expansion > ExpansionSharedPtr
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
boost::shared_ptr< GeometryVector > Composite
bool IsLeftAdjacentEdge(const int n, const int e)
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
static const NekDouble kNekUnsetDouble
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.
AssemblyMapDGSharedPtr m_traceMap
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
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.
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::CommSharedPtr m_comm
Communicator.
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
GlobalLinSysMapShPtr m_globalBndMat
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
virtual void v_GetFwdBwdTracePhys(const Array< OneD, const NekDouble > &field, 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...
boost::shared_ptr< TriGeom > TriGeomSharedPtr
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
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.
const BoundaryRegionCollection & GetBoundaryRegions(void) const
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...
Array< OneD, int > m_BCtoEdgMap
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
boost::shared_ptr< Basis > BasisSharedPtr
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
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph2D, const SpatialDomains::BoundaryConditions &bcs, const std::string &variable, const bool DeclareCoeffPhysArrays=true)
This function discretises the boundary conditions by setting up a list of one-dimensional boundary ex...
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Set up a list of element IDs and edge IDs that link to the boundary conditions.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
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.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr