Skip to content

GitLab

Commit f06db8fb authored by Douglas Serson's avatar Douglas Serson
Browse files

Merge branch 'fix/GetOffsetElmtId' into 'master' 

Fix/rm_getoffsetelmtid

See merge request !758
parents 9738e4ae 4b560009
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
View file @ f06db8fb
......@@ -17,6 +17,8 @@ v4.5.0
v4.4.1
------
**Library**
- Remove m_offset_elmt_id and GetOffsetElmtId which fixed problems in 2D when
quad elements are listed before tri elements (!758)
- Remove the duplicate output of errorutil (!756)
- Fix interpolation issue with Lagrange basis functions (!768)
......
library/Demos/MultiRegions/Tests/Helmholtz3D_CG_Prism_Deformed.tst
View file @ f06db8fb
......@@ -9,10 +9,10 @@
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-9">3.35253e-06</value>
<value tolerance="1e-9">1.13842e-05</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-9">4.45606e-05</value>
<value tolerance="1e-9">7.59576e-05</value>
</metric>
</metrics>
</test>
......
library/FieldUtils/OutputModules/OutputTecplot.cpp
View file @ f06db8fb
......@@ -219,7 +219,6 @@ void OutputTecplot::Process(po::variables_map &vm)
MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
m_coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
var.insert(var.begin(), coordVars, coordVars + m_coordim);
if (HomoExpType == MultiRegions::e3DH1D)
{
......@@ -242,6 +241,8 @@ void OutputTecplot::Process(po::variables_map &vm)
m_coordim += 2;
}
var.insert(var.begin(), coordVars, coordVars + m_coordim);
m_zoneType = (TecplotZoneType)(2*(nBases-1) + 1);
// Calculate connectivity
......
library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
View file @ f06db8fb
......@@ -275,7 +275,7 @@ namespace Nektar
for(j = 0; j < locExpVector.size(); j++)
{
exp = locExpVector[locExp.GetOffset_Elmt_Id(j)];
exp = locExpVector[j];
for(k = 0; k < exp->GetNverts(); k++)
{
......@@ -580,7 +580,7 @@ namespace Nektar
/// - Count verts, edges, face and add up edges and face sizes
for(i = 0; i < locExpVector.size(); ++i)
{
exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
exp = locExpVector[i];
nTotalVerts += exp->GetNverts();
nTotalEdges += exp->GetNedges();
nTotalFaces += exp->GetNfaces();
......@@ -1073,7 +1073,7 @@ namespace Nektar
// edges respectively which are local to this process.
for(i = cnt = 0; i < locExpVector.size(); ++i)
{
int elmtid = locExp.GetOffset_Elmt_Id(i);
int elmtid = i;
exp = locExpVector[elmtid];
for (j = 0; j < exp->GetNverts(); ++j)
{
......@@ -1552,7 +1552,7 @@ namespace Nektar
for(i = 0; i < locExpVector.size(); ++i)
{
exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
exp = locExpVector[i];
for(j = 0; j < exp->GetNverts(); ++j)
{
......@@ -1612,10 +1612,10 @@ namespace Nektar
for(i = 0; i < m_numPatches; ++i)
{
m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
locExpVector[i]->NumBndryCoeffs();
m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
locExpVector[i]->GetNcoeffs() -
locExpVector[i]->NumBndryCoeffs();
}
/**
......@@ -1987,7 +1987,7 @@ namespace Nektar
for (i = 0; i < locExpVector.size(); ++i)
{
exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
exp = locExpVector[i];
for (j = 0; j < exp->GetNverts(); ++j)
{
......
library/MultiRegions/AssemblyMap/AssemblyMapDG.cpp
View file @ f06db8fb
......@@ -178,7 +178,7 @@ namespace Nektar
int nbndry = 0;
for(i = 0; i < nel; ++i) // count number of elements in array
{
eid = locExp.GetOffset_Elmt_Id(i);
eid = i;
nbndry += expList[eid]->NumDGBndryCoeffs();
m_numLocalIntCoeffsPerPatch[i] = 0;
m_numLocalBndCoeffsPerPatch[i] =
......@@ -230,7 +230,7 @@ namespace Nektar
// Set up boost Graph
for(i = 0; i < nel; ++i)
{
eid = locExp.GetOffset_Elmt_Id(i);
eid = i;
for(j = 0; j < expList[eid]->GetNtrace(); ++j)
{
......@@ -319,8 +319,7 @@ namespace Nektar
cnt = 0;
for(i = 0; i < nel; ++i)
{
// order list according to m_offset_elmt_id details in expList
eid = locExp.GetOffset_Elmt_Id(i);
eid = i;
exp = expList[eid];
for(j = 0; j < exp->GetNtrace(); ++j)
......@@ -585,9 +584,7 @@ namespace Nektar
cnt = 0;
for(i = 0; i < locExpVector.size(); ++i)
{
// Order list according to m_offset_elmt_id details in Exp
// so that triangules are listed first and then quads
eid = locExp.GetOffset_Elmt_Id(i);
eid = i;
locExpansion = locExpVector[eid];
nDim = locExpansion->GetShapeDimension();
......
library/MultiRegions/DisContField2D.cpp
View file @ f06db8fb
......@@ -1814,7 +1814,7 @@ namespace Nektar
StdRegions::Orientation edgedir;
int eid,cnt;
int cnt;
int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
......@@ -1824,8 +1824,6 @@ namespace Nektar
// Calculate Q using standard DG formulation.
for(i = cnt = 0; i < GetExpSize(); ++i)
{
eid = m_offset_elmt_id[i];
// Probably a better way of setting up lambda than this.
// Note cannot use PutCoeffsInToElmts since lambda space
// is mapped during the solve.
......@@ -1834,21 +1832,21 @@ namespace Nektar
for(e = 0; e < nEdges; ++e)
{
edgedir = (*m_exp)[eid]->GetEorient(e);
ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs();
edgedir = (*m_exp)[i]->GetEorient(e);
ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
elmtToTrace[eid][e]->SetCoeffsToOrientation(
elmtToTrace[i][e]->SetCoeffsToOrientation(
edgedir, edgeCoeffs[e], edgeCoeffs[e]);
edge_lambda = edge_lambda + ncoeff_edge;
}
(*m_exp)[eid]->DGDeriv(dir,
tmp_coeffs=m_coeffs+m_coeff_offset[eid],
elmtToTrace[eid],
(*m_exp)[i]->DGDeriv(dir,
tmp_coeffs=m_coeffs+m_coeff_offset[i],
elmtToTrace[i],
edgeCoeffs,
out_tmp = out_d+cnt);
cnt += (*m_exp)[eid]->GetNcoeffs();
cnt += (*m_exp)[i]->GetNcoeffs();
}
BwdTrans(out_d,m_phys);
......@@ -1923,9 +1921,9 @@ namespace Nektar
// Loop over all expansions in the domain
for(cnt = cnt1 = n = 0; n < nexp; ++n)
{
nbndry = (*m_exp)[m_offset_elmt_id[n]]->NumDGBndryCoeffs();
nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
e_ncoeffs = (*m_exp)[m_offset_elmt_id[n]]->GetNcoeffs();
e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
e_f = f + cnt;
e_l = loc_lambda + cnt1;
......@@ -2123,7 +2121,7 @@ namespace Nektar
StdRegions::Orientation edgedir;
int eid,nq_elmt, nm_elmt;
int nq_elmt, nm_elmt;
int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
Array<OneD, NekDouble> tmp_coeffs;
......@@ -2134,20 +2132,18 @@ namespace Nektar
// Calculate Q using standard DG formulation.
for(i = cnt = 0; i < GetExpSize(); ++i)
{
eid = m_offset_elmt_id[i];
nq_elmt = (*m_exp)[eid]->GetTotPoints();
nm_elmt = (*m_exp)[eid]->GetNcoeffs();
nq_elmt = (*m_exp)[i]->GetTotPoints();
nm_elmt = (*m_exp)[i]->GetNcoeffs();
qrhs = Array<OneD, NekDouble>(nq_elmt);
qrhs1 = Array<OneD, NekDouble>(nq_elmt);
force = Array<OneD, NekDouble>(2*nm_elmt);
out_tmp = force + nm_elmt;
LocalRegions::ExpansionSharedPtr ppExp;
int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints();
int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints();
int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes();
int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes();
int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
// Probably a better way of setting up lambda than this. Note
// cannot use PutCoeffsInToElmts since lambda space is mapped
......@@ -2155,19 +2151,19 @@ namespace Nektar
int nEdges = (*m_exp)[i]->GetNedges();
Array<OneD, Array<OneD, NekDouble> > edgeCoeffs(nEdges);
for(e = 0; e < (*m_exp)[eid]->GetNedges(); ++e)
for(e = 0; e < (*m_exp)[i]->GetNedges(); ++e)
{
edgedir = (*m_exp)[eid]->GetEorient(e);
ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs();
edgedir = (*m_exp)[i]->GetEorient(e);
ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
elmtToTrace[eid][e]->SetCoeffsToOrientation(
elmtToTrace[i][e]->SetCoeffsToOrientation(
edgedir, edgeCoeffs[e], edgeCoeffs[e]);
edge_lambda = edge_lambda + ncoeff_edge;
}
//creating orthogonal expansion (checking if we have quads or triangles)
LibUtilities::ShapeType shape = (*m_exp)[eid]->DetShapeType();
LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
switch(shape)
{
case LibUtilities::eQuadrilateral:
......@@ -2176,7 +2172,7 @@ namespace Nektar
const LibUtilities::PointsKey PkeyQ2(num_points1,LibUtilities::eGaussLobattoLegendre);
LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A, num_modes0, PkeyQ1);
LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A, num_modes1, PkeyQ2);
SpatialDomains::QuadGeomSharedPtr qGeom = boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>((*m_exp)[eid]->GetGeom());
SpatialDomains::QuadGeomSharedPtr qGeom = boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>((*m_exp)[i]->GetGeom());
ppExp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(BkeyQ1, BkeyQ2, qGeom);
}
break;
......@@ -2186,7 +2182,7 @@ namespace Nektar
const LibUtilities::PointsKey PkeyT2(num_points1,LibUtilities::eGaussRadauMAlpha1Beta0);
LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1);
LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes1, PkeyT2);
SpatialDomains::TriGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TriGeom>((*m_exp)[eid]->GetGeom());
SpatialDomains::TriGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TriGeom>((*m_exp)[i]->GetGeom());
ppExp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(BkeyT1, BkeyT2, tGeom);
}
break;
......@@ -2197,29 +2193,29 @@ namespace Nektar
//DGDeriv
// (d/dx w, d/dx q_0)
(*m_exp)[eid]->DGDeriv(
0,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
elmtToTrace[eid], edgeCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[eid]->IProductWRTDerivBase(0,qrhs,force);
(*m_exp)[i]->DGDeriv(
0,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[i], edgeCoeffs, out_tmp);
(*m_exp)[i]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[i]->IProductWRTDerivBase(0,qrhs,force);
ppExp->IProductWRTDerivBase(0,qrhs,force);
// + (d/dy w, d/dy q_1)
(*m_exp)[eid]->DGDeriv(
1,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
elmtToTrace[eid], edgeCoeffs, out_tmp);
(*m_exp)[i]->DGDeriv(
1,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[i], edgeCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[eid]->IProductWRTDerivBase(1,qrhs,out_tmp);
(*m_exp)[i]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[i]->IProductWRTDerivBase(1,qrhs,out_tmp);
ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
// determine force[0] = (1,u)
(*m_exp)[eid]->BwdTrans(
tmp_coeffs = m_coeffs + m_coeff_offset[eid],qrhs);
force[0] = (*m_exp)[eid]->Integral(qrhs);
(*m_exp)[i]->BwdTrans(
tmp_coeffs = m_coeffs + m_coeff_offset[i],qrhs);
force[0] = (*m_exp)[i]->Integral(qrhs);
// multiply by inverse Laplacian matrix
// get matrix inverse
......@@ -2234,7 +2230,7 @@ namespace Nektar
// Transforming back to modified basis
Array<OneD, NekDouble> work(nq_elmt);
ppExp->BwdTrans(out.GetPtr(), work);
(*m_exp)[eid]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[eid]);
(*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
}
}
......
library/MultiRegions/DisContField3D.cpp
View file @ f06db8fb
......@@ -2081,9 +2081,9 @@ using namespace boost::assign;
// Loop over all expansions in the domain
for(cnt = cnt1 = n = 0; n < nexp; ++n)
{
nbndry = (*m_exp)[m_offset_elmt_id[n]]->NumDGBndryCoeffs();
nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
e_ncoeffs = (*m_exp)[m_offset_elmt_id[n]]->GetNcoeffs();
e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
e_f = f + cnt;
e_l = loc_lambda + cnt1;
......@@ -2273,8 +2273,8 @@ using namespace boost::assign;
Array<OneD, NekDouble> force, out_tmp,qrhs,qrhs1;
Array<OneD, Array< OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
int eid,nq_elmt, nm_elmt;
int nq_elmt, nm_elmt;
int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), face_lambda;
Array<OneD, NekDouble> tmp_coeffs;
......@@ -2288,30 +2288,29 @@ using namespace boost::assign;
LocalRegions::Expansion3DSharedPtr exp =
(*m_exp)[i]->as<LocalRegions::Expansion3D>();
eid = m_offset_elmt_id[i];
nq_elmt = (*m_exp)[eid]->GetTotPoints();
nm_elmt = (*m_exp)[eid]->GetNcoeffs();
nq_elmt = (*m_exp)[i]->GetTotPoints();
nm_elmt = (*m_exp)[i]->GetNcoeffs();
qrhs = Array<OneD, NekDouble>(nq_elmt);
qrhs1 = Array<OneD, NekDouble>(nq_elmt);
force = Array<OneD, NekDouble>(2*nm_elmt);
out_tmp = force + nm_elmt;
LocalRegions::ExpansionSharedPtr ppExp;
int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints();
int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints();
int num_points2 = (*m_exp)[eid]->GetBasis(2)->GetNumPoints();
int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes();
int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes();
int num_modes2 = (*m_exp)[eid]->GetBasis(2)->GetNumModes();
int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
int num_points2 = (*m_exp)[i]->GetBasis(2)->GetNumPoints();
int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
int num_modes2 = (*m_exp)[i]->GetBasis(2)->GetNumModes();
// Probably a better way of setting up lambda than this. Note
// cannot use PutCoeffsInToElmts since lambda space is mapped
// during the solve.
int nFaces = (*m_exp)[eid]->GetNfaces();
int nFaces = (*m_exp)[i]->GetNfaces();
Array<OneD, Array<OneD, NekDouble> > faceCoeffs(nFaces);
for(f = 0; f < nFaces; ++f)
{
ncoeff_face = elmtToTrace[eid][f]->GetNcoeffs();
ncoeff_face = elmtToTrace[i][f]->GetNcoeffs();
faceCoeffs[f] = Array<OneD, NekDouble>(ncoeff_face);
Vmath::Vcopy(ncoeff_face, face_lambda, 1, faceCoeffs[f], 1);
exp->SetFaceToGeomOrientation(f, faceCoeffs[f]);
......@@ -2319,7 +2318,7 @@ using namespace boost::assign;
}
//creating orthogonal expansion (checking if we have quads or triangles)
LibUtilities::ShapeType shape = (*m_exp)[eid]->DetShapeType();
LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
switch(shape)
{
case LibUtilities::eHexahedron:
......@@ -2330,7 +2329,7 @@ using namespace boost::assign;
LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A, num_modes0, PkeyH1);
LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A, num_modes1, PkeyH2);
LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A, num_modes2, PkeyH3);
SpatialDomains::HexGeomSharedPtr hGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>((*m_exp)[eid]->GetGeom());
SpatialDomains::HexGeomSharedPtr hGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>((*m_exp)[i]->GetGeom());
ppExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(BkeyH1, BkeyH2, BkeyH3, hGeom);
}
break;
......@@ -2342,7 +2341,7 @@ using namespace boost::assign;
LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1);
LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes1, PkeyT2);
LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C, num_modes2, PkeyT3);
SpatialDomains::TetGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>((*m_exp)[eid]->GetGeom());
SpatialDomains::TetGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>((*m_exp)[i]->GetGeom());
ppExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(BkeyT1, BkeyT2, BkeyT3, tGeom);
}
break;
......@@ -2354,7 +2353,7 @@ using namespace boost::assign;
LibUtilities::BasisKey BkeyP1(LibUtilities::eOrtho_A, num_modes0, PkeyP1);
LibUtilities::BasisKey BkeyP2(LibUtilities::eOrtho_A, num_modes1, PkeyP2);
LibUtilities::BasisKey BkeyP3(LibUtilities::eOrtho_B, num_modes2, PkeyP3);
SpatialDomains::PrismGeomSharedPtr pGeom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>((*m_exp)[eid]->GetGeom());
SpatialDomains::PrismGeomSharedPtr pGeom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>((*m_exp)[i]->GetGeom());
ppExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(BkeyP1, BkeyP2, BkeyP3, pGeom);
}
break;
......@@ -2365,34 +2364,34 @@ using namespace boost::assign;
//DGDeriv
// (d/dx w, q_0)
(*m_exp)[eid]->DGDeriv(
0,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
elmtToTrace[eid], faceCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
(*m_exp)[i]->DGDeriv(
0,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[i], faceCoeffs, out_tmp);
(*m_exp)[i]->BwdTrans(out_tmp,qrhs);
ppExp->IProductWRTDerivBase(0,qrhs,force);
// + (d/dy w, q_1)
(*m_exp)[eid]->DGDeriv(
1,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
elmtToTrace[eid], faceCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
(*m_exp)[i]->DGDeriv(
1,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[i], faceCoeffs, out_tmp);
(*m_exp)[i]->BwdTrans(out_tmp,qrhs);
ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
// + (d/dz w, q_2)
(*m_exp)[eid]->DGDeriv(
2,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
elmtToTrace[eid], faceCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
(*m_exp)[i]->DGDeriv(
2,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[i], faceCoeffs, out_tmp);
(*m_exp)[i]->BwdTrans(out_tmp,qrhs);
ppExp->IProductWRTDerivBase(2,qrhs,out_tmp);
Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
// determine force[0] = (1,u)
(*m_exp)[eid]->BwdTrans(
tmp_coeffs = m_coeffs + m_coeff_offset[eid],qrhs);
force[0] = (*m_exp)[eid]->Integral(qrhs);
(*m_exp)[i]->BwdTrans(
tmp_coeffs = m_coeffs + m_coeff_offset[i],qrhs);
force[0] = (*m_exp)[i]->Integral(qrhs);
// multiply by inverse Laplacian matrix
// get matrix inverse
......@@ -2407,8 +2406,8 @@ using namespace boost::assign;
// Transforming back to modified basis
Array<OneD, NekDouble> work(nq_elmt);
ppExp->BwdTrans(out.GetPtr(), work);
(*m_exp)[eid]->FwdTrans(work,
tmp_coeffs = outarray + m_coeff_offset[eid]);
(*m_exp)[i]->FwdTrans(work,
tmp_coeffs = outarray + m_coeff_offset[i]);
}
}
......
library/MultiRegions/ExpList.cpp
View file @ f06db8fb
......@@ -105,7 +105,6 @@ namespace Nektar
::AllocateSharedPtr()),
m_coeff_offset(),
m_phys_offset(),
m_offset_elmt_id(),
m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
m_WaveSpace(false)
{
......@@ -132,7 +131,6 @@ namespace Nektar
::AllocateSharedPtr()),
m_coeff_offset(),
m_phys_offset(),
m_offset_elmt_id(),
m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
m_WaveSpace(false)
{
......@@ -160,7 +158,6 @@ namespace Nektar
::AllocateSharedPtr()),
m_coeff_offset(),
m_phys_offset(),
m_offset_elmt_id(),
m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
m_WaveSpace(false)
{
......@@ -187,7 +184,6 @@ namespace Nektar
::AllocateSharedPtr()),
m_coeff_offset(),
m_phys_offset(),
m_offset_elmt_id(),
m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
m_WaveSpace(false)
{
......@@ -225,7 +221,6 @@ namespace Nektar
m_coll_phys_offset(in.m_coll_phys_offset),
m_coeff_offset(in.m_coeff_offset),
m_phys_offset(in.m_phys_offset),
m_offset_elmt_id(in.m_offset_elmt_id),
m_globalOptParam(in.m_globalOptParam),
m_blockMat(in.m_blockMat),
m_WaveSpace(false)
......@@ -239,6 +234,35 @@ namespace Nektar
}
}
/**
* Each expansion (local element) is processed in turn to
* determine the number of coefficients and physical data
* points it contributes to the domain. Twoe arrays,
* #m_coeff_offset are #m_phys_offset are also initialised and
* updated to store the data offsets of each element in the
* #m_coeffs and #m_phys arrays, and the element id that each
* consecutive block is associated respectively.
*/
void ExpList::SetCoeffPhysOffsets()
{
int i;
// Set up offset information and array sizes
m_coeff_offset = Array<OneD,int>(m_exp->size());
m_phys_offset = Array<OneD,int>(m_exp->size());
m_ncoeffs = m_npoints = 0;
for(i = 0; i < m_exp->size(); ++i)
{
m_coeff_offset[i] = m_ncoeffs;
m_phys_offset [i] = m_npoints;
m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
m_npoints += (*m_exp)[i]->GetTotPoints();
}
}
/**
*
*/
......@@ -765,10 +789,10 @@ namespace Nektar
{
for(i = 0 ; i < (*m_exp).size(); ++i)
{
if((*m_exp)[m_offset_elmt_id[i]]->DetShapeType()
if((*m_exp)[i]->DetShapeType()
== ShapeType)
{
elmt_id[n_exp++] = m_offset_elmt_id[i];
elmt_id[n_exp++] = i;
}
}
}
......@@ -777,7 +801,7 @@ namespace Nektar
n_exp = (*m_exp).size();
for(i = 0; i < n_exp; ++i)
{
elmt_id[i] = m_offset_elmt_id[i];
elmt_id[i] = i;
}
}
......@@ -922,9 +946,9 @@ namespace Nektar