Inconsistent strategy in dealing with 1D and 2D/3D expansions in DisContField::v_GetBoundaryToElmtMap
Here is a snap of DisContField.cpp:
// Set up a list of element ids and trace ids that link to the
// boundary conditions
void DisContField::v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
Array<OneD, int> &TraceID)
{
if (m_BCtoElmMap.size() == 0)
{
switch (m_expType)
{
case e1D:
{
map<int, int> VertGID;
int i, n, id;
int bid, cnt, Vid;
int nbcs = m_bndConditions.size();
// make sure arrays are of sufficient length
m_BCtoElmMap = Array<OneD, int>(nbcs, -1);
m_BCtoTraceMap = Array<OneD, int>(nbcs);
...
}
case e2D:
{
map<int, int> globalIdMap;
int i, n;
int cnt;
int nbcs = 0;
// Populate global ID map (takes global geometry
// ID to local expansion list ID).
for (i = 0; i < GetExpSize(); ++i)
{
globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
}
// Determine number of boundary condition expansions.
for (i = 0; i < m_bndConditions.size(); ++i)
{
nbcs += m_bndCondExpansions[i]->GetExpSize();
}
// Initialize arrays
m_BCtoElmMap = Array<OneD, int>(nbcs);
m_BCtoTraceMap = Array<OneD, int>(nbcs);
...
For e1D case, the size of array m_BCtoElmMap
and m_BCtoTraceMap
is m_bndConditions.size
. However, for e2D and e3D cases, the size becomes sum of all bndcondition expansions.
In 1D problem, the boundary regions are typically the 2 end points of a line mesh. For example
<COMPOSITE>
<C ID="0"> S[0-9] </C>
<C ID="1"> V[0] </C>
<C ID="2"> V[10] </C>
</COMPOSITE>
If we try to include 2 endpoints in single boundary region, (see solvers/AcousticSolver/Tests/APE_1DPulseSource_WeakDG_MODIFIED.xml)
<BOUNDARYREGIONS>
<B ID="0"> C[1,2] </B>
</BOUNDARYREGIONS>
Then the following part of the code will give error because the m_bndCondExpansions[i]->locExpList
contains more than 1 expansion thus the subscript [cnt + e]
will exceed the bound of array TraceID
.
if (SetUpJustDG)
{
SetUpDG(variable, ImpType);
}
else
{
int i, cnt;
Array<OneD, int> ElmtID, TraceID;
GetBoundaryToElmtMap(ElmtID, TraceID);
for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
{
MultiRegions::ExpListSharedPtr locExpList;
locExpList = m_bndCondExpansions[i];
for (int e = 0; e < locExpList->GetExpSize(); ++e)
{
(*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
locExpList->GetExp(e));
locExpList->GetExp(e)->SetAdjacentElementExp(
TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
}
cnt += m_bndCondExpansions[i]->GetExpSize();
}
But, if we define 2 boundary regions and each contains only 1 end point, namely:
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
Then the above code works fine and gives no error. But it is just a coincidence I have to say.
At present, nearly all of the 1D test cases follow the one-point-in-one-region rule, except APE_1DPulseSource_WeakDG_MODIFIED.
But luckily, that case enables SetUpJustDG
and the above part of code is not executed! That is why we do not receive any error so far!
Recently, I am trying to add this part of code to the SetUpJustDG
branch for other features, and just happened to detect this issue.
I post it in the hope that everyone can know it and make consensus on how to solve it.
I am new to Nektar group and thanks to all.