Commit 7fb9f05d authored by Spencer Sherwin's avatar Spencer Sherwin

Removed "eid = i" statements from previous removal of GetOffset_ElmtId calls

parent 24d28922
...@@ -1814,7 +1814,7 @@ namespace Nektar ...@@ -1814,7 +1814,7 @@ namespace Nektar
StdRegions::Orientation edgedir; StdRegions::Orientation edgedir;
int eid,cnt; int cnt;
int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs(); int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda; Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda); m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
...@@ -1824,8 +1824,6 @@ namespace Nektar ...@@ -1824,8 +1824,6 @@ namespace Nektar
// Calculate Q using standard DG formulation. // Calculate Q using standard DG formulation.
for(i = cnt = 0; i < GetExpSize(); ++i) for(i = cnt = 0; i < GetExpSize(); ++i)
{ {
eid = i;
// Probably a better way of setting up lambda than this. // Probably a better way of setting up lambda than this.
// Note cannot use PutCoeffsInToElmts since lambda space // Note cannot use PutCoeffsInToElmts since lambda space
// is mapped during the solve. // is mapped during the solve.
...@@ -1834,21 +1832,21 @@ namespace Nektar ...@@ -1834,21 +1832,21 @@ namespace Nektar
for(e = 0; e < nEdges; ++e) for(e = 0; e < nEdges; ++e)
{ {
edgedir = (*m_exp)[eid]->GetEorient(e); edgedir = (*m_exp)[i]->GetEorient(e);
ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs(); ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge); edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1); Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
elmtToTrace[eid][e]->SetCoeffsToOrientation( elmtToTrace[i][e]->SetCoeffsToOrientation(
edgedir, edgeCoeffs[e], edgeCoeffs[e]); edgedir, edgeCoeffs[e], edgeCoeffs[e]);
edge_lambda = edge_lambda + ncoeff_edge; edge_lambda = edge_lambda + ncoeff_edge;
} }
(*m_exp)[eid]->DGDeriv(dir, (*m_exp)[i]->DGDeriv(dir,
tmp_coeffs=m_coeffs+m_coeff_offset[eid], tmp_coeffs=m_coeffs+m_coeff_offset[i],
elmtToTrace[eid], elmtToTrace[i],
edgeCoeffs, edgeCoeffs,
out_tmp = out_d+cnt); out_tmp = out_d+cnt);
cnt += (*m_exp)[eid]->GetNcoeffs(); cnt += (*m_exp)[i]->GetNcoeffs();
} }
BwdTrans(out_d,m_phys); BwdTrans(out_d,m_phys);
...@@ -2123,7 +2121,7 @@ namespace Nektar ...@@ -2123,7 +2121,7 @@ namespace Nektar
StdRegions::Orientation edgedir; StdRegions::Orientation edgedir;
int eid,nq_elmt, nm_elmt; int nq_elmt, nm_elmt;
int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs(); int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda; Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
Array<OneD, NekDouble> tmp_coeffs; Array<OneD, NekDouble> tmp_coeffs;
...@@ -2134,20 +2132,18 @@ namespace Nektar ...@@ -2134,20 +2132,18 @@ namespace Nektar
// Calculate Q using standard DG formulation. // Calculate Q using standard DG formulation.
for(i = cnt = 0; i < GetExpSize(); ++i) for(i = cnt = 0; i < GetExpSize(); ++i)
{ {
eid = i; nq_elmt = (*m_exp)[i]->GetTotPoints();
nm_elmt = (*m_exp)[i]->GetNcoeffs();
nq_elmt = (*m_exp)[eid]->GetTotPoints();
nm_elmt = (*m_exp)[eid]->GetNcoeffs();
qrhs = Array<OneD, NekDouble>(nq_elmt); qrhs = Array<OneD, NekDouble>(nq_elmt);
qrhs1 = Array<OneD, NekDouble>(nq_elmt); qrhs1 = Array<OneD, NekDouble>(nq_elmt);
force = Array<OneD, NekDouble>(2*nm_elmt); force = Array<OneD, NekDouble>(2*nm_elmt);
out_tmp = force + nm_elmt; out_tmp = force + nm_elmt;
LocalRegions::ExpansionSharedPtr ppExp; LocalRegions::ExpansionSharedPtr ppExp;
int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints(); int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints(); int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes(); int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes(); int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
// Probably a better way of setting up lambda than this. Note // Probably a better way of setting up lambda than this. Note
// cannot use PutCoeffsInToElmts since lambda space is mapped // cannot use PutCoeffsInToElmts since lambda space is mapped
...@@ -2155,19 +2151,19 @@ namespace Nektar ...@@ -2155,19 +2151,19 @@ namespace Nektar
int nEdges = (*m_exp)[i]->GetNedges(); int nEdges = (*m_exp)[i]->GetNedges();
Array<OneD, Array<OneD, NekDouble> > edgeCoeffs(nEdges); 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); edgedir = (*m_exp)[i]->GetEorient(e);
ncoeff_edge = elmtToTrace[eid][e]->GetNcoeffs(); ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge); edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1); Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
elmtToTrace[eid][e]->SetCoeffsToOrientation( elmtToTrace[i][e]->SetCoeffsToOrientation(
edgedir, edgeCoeffs[e], edgeCoeffs[e]); edgedir, edgeCoeffs[e], edgeCoeffs[e]);
edge_lambda = edge_lambda + ncoeff_edge; edge_lambda = edge_lambda + ncoeff_edge;
} }
//creating orthogonal expansion (checking if we have quads or triangles) //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) switch(shape)
{ {
case LibUtilities::eQuadrilateral: case LibUtilities::eQuadrilateral:
...@@ -2176,7 +2172,7 @@ namespace Nektar ...@@ -2176,7 +2172,7 @@ namespace Nektar
const LibUtilities::PointsKey PkeyQ2(num_points1,LibUtilities::eGaussLobattoLegendre); const LibUtilities::PointsKey PkeyQ2(num_points1,LibUtilities::eGaussLobattoLegendre);
LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A, num_modes0, PkeyQ1); LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A, num_modes0, PkeyQ1);
LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A, num_modes1, PkeyQ2); 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); ppExp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(BkeyQ1, BkeyQ2, qGeom);
} }
break; break;
...@@ -2186,7 +2182,7 @@ namespace Nektar ...@@ -2186,7 +2182,7 @@ namespace Nektar
const LibUtilities::PointsKey PkeyT2(num_points1,LibUtilities::eGaussRadauMAlpha1Beta0); const LibUtilities::PointsKey PkeyT2(num_points1,LibUtilities::eGaussRadauMAlpha1Beta0);
LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1); LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1);
LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes1, PkeyT2); 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); ppExp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(BkeyT1, BkeyT2, tGeom);
} }
break; break;
...@@ -2197,29 +2193,29 @@ namespace Nektar ...@@ -2197,29 +2193,29 @@ namespace Nektar
//DGDeriv //DGDeriv
// (d/dx w, d/dx q_0) // (d/dx w, d/dx q_0)
(*m_exp)[eid]->DGDeriv( (*m_exp)[i]->DGDeriv(
0,tmp_coeffs = m_coeffs + m_coeff_offset[eid], 0,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[eid], edgeCoeffs, out_tmp); elmtToTrace[i], edgeCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs); (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[eid]->IProductWRTDerivBase(0,qrhs,force); //(*m_exp)[i]->IProductWRTDerivBase(0,qrhs,force);
ppExp->IProductWRTDerivBase(0,qrhs,force); ppExp->IProductWRTDerivBase(0,qrhs,force);
// + (d/dy w, d/dy q_1) // + (d/dy w, d/dy q_1)
(*m_exp)[eid]->DGDeriv( (*m_exp)[i]->DGDeriv(
1,tmp_coeffs = m_coeffs + m_coeff_offset[eid], 1,tmp_coeffs = m_coeffs + m_coeff_offset[i],
elmtToTrace[eid], edgeCoeffs, out_tmp); elmtToTrace[i], edgeCoeffs, out_tmp);
(*m_exp)[eid]->BwdTrans(out_tmp,qrhs); (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
//(*m_exp)[eid]->IProductWRTDerivBase(1,qrhs,out_tmp); //(*m_exp)[i]->IProductWRTDerivBase(1,qrhs,out_tmp);
ppExp->IProductWRTDerivBase(1,qrhs,out_tmp); ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1); Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
// determine force[0] = (1,u) // determine force[0] = (1,u)
(*m_exp)[eid]->BwdTrans( (*m_exp)[i]->BwdTrans(
tmp_coeffs = m_coeffs + m_coeff_offset[eid],qrhs); tmp_coeffs = m_coeffs + m_coeff_offset[i],qrhs);
force[0] = (*m_exp)[eid]->Integral(qrhs); force[0] = (*m_exp)[i]->Integral(qrhs);
// multiply by inverse Laplacian matrix // multiply by inverse Laplacian matrix
// get matrix inverse // get matrix inverse
...@@ -2234,7 +2230,7 @@ namespace Nektar ...@@ -2234,7 +2230,7 @@ namespace Nektar
// Transforming back to modified basis // Transforming back to modified basis
Array<OneD, NekDouble> work(nq_elmt); Array<OneD, NekDouble> work(nq_elmt);
ppExp->BwdTrans(out.GetPtr(), work); 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]);
} }
} }
......
...@@ -2359,18 +2359,17 @@ namespace Nektar ...@@ -2359,18 +2359,17 @@ namespace Nektar
{ {
std::vector<unsigned int> nummodes; std::vector<unsigned int> nummodes;
vector<LibUtilities::BasisType> basisTypes; vector<LibUtilities::BasisType> basisTypes;
int eid = i; for(int j= 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
{ {
nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j)); nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
basisTypes.push_back(fromExpList->GetExp(eid)->GetBasisType(j)); basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
} }
(*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0, (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
&toCoeffs[m_coeff_offset[eid]], &toCoeffs[m_coeff_offset[i]],
basisTypes); basisTypes);
offset += fromExpList->GetExp(eid)->GetNcoeffs(); offset += fromExpList->GetExp(i)->GetNcoeffs();
} }
} }
......
...@@ -1227,20 +1227,19 @@ namespace Nektar ...@@ -1227,20 +1227,19 @@ namespace Nektar
for(int i = 0; i < GetExpSize(); ++i) for(int i = 0; i < GetExpSize(); ++i)
{ {
// get new points key // get new points key
int eid = i; int pt0 = (*m_exp)[i]->GetNumPoints(0);
int pt0 = (*m_exp)[eid]->GetNumPoints(0); int pt1 = (*m_exp)[i]->GetNumPoints(1);
int pt1 = (*m_exp)[eid]->GetNumPoints(1);
int npt0 = (int) pt0*scale; int npt0 = (int) pt0*scale;
int npt1 = (int) pt1*scale; int npt1 = (int) pt1*scale;
LibUtilities::PointsKey newPointsKey0(npt0, LibUtilities::PointsKey newPointsKey0(npt0,
(*m_exp)[eid]->GetPointsType(0)); (*m_exp)[i]->GetPointsType(0));
LibUtilities::PointsKey newPointsKey1(npt1, LibUtilities::PointsKey newPointsKey1(npt1,
(*m_exp)[eid]->GetPointsType(1)); (*m_exp)[i]->GetPointsType(1));
// Interpolate points; // Interpolate points;
LibUtilities::Interp2D((*m_exp)[eid]->GetBasis(0)->GetPointsKey(), LibUtilities::Interp2D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(1)->GetPointsKey(), (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
&inarray[cnt],newPointsKey0, &inarray[cnt],newPointsKey0,
newPointsKey1,&outarray[cnt1]); newPointsKey1,&outarray[cnt1]);
...@@ -1260,24 +1259,23 @@ namespace Nektar ...@@ -1260,24 +1259,23 @@ namespace Nektar
for(int i = 0; i < GetExpSize(); ++i) for(int i = 0; i < GetExpSize(); ++i)
{ {
// get new points key // get new points key
int eid = i; int pt0 = (*m_exp)[i]->GetNumPoints(0);
int pt0 = (*m_exp)[eid]->GetNumPoints(0); int pt1 = (*m_exp)[i]->GetNumPoints(1);
int pt1 = (*m_exp)[eid]->GetNumPoints(1);
int npt0 = (int) pt0*scale; int npt0 = (int) pt0*scale;
int npt1 = (int) pt1*scale; int npt1 = (int) pt1*scale;
LibUtilities::PointsKey newPointsKey0(npt0, LibUtilities::PointsKey newPointsKey0(npt0,
(*m_exp)[eid]->GetPointsType(0)); (*m_exp)[i]->GetPointsType(0));
LibUtilities::PointsKey newPointsKey1(npt1, LibUtilities::PointsKey newPointsKey1(npt1,
(*m_exp)[eid]->GetPointsType(1)); (*m_exp)[i]->GetPointsType(1));
// Project points; // Project points;
LibUtilities::PhysGalerkinProject2D( LibUtilities::PhysGalerkinProject2D(
newPointsKey0, newPointsKey0,
newPointsKey1, newPointsKey1,
&inarray[cnt], &inarray[cnt],
(*m_exp)[eid]->GetBasis(0)->GetPointsKey(), (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(1)->GetPointsKey(), (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
&outarray[cnt1]); &outarray[cnt1]);
cnt += npt0*npt1; cnt += npt0*npt1;
......
...@@ -572,22 +572,21 @@ namespace Nektar ...@@ -572,22 +572,21 @@ namespace Nektar
for(int i = 0; i < GetExpSize(); ++i) for(int i = 0; i < GetExpSize(); ++i)
{ {
// get new points key // get new points key
int eid = i; int pt0 = (*m_exp)[i]->GetNumPoints(0);
int pt0 = (*m_exp)[eid]->GetNumPoints(0); int pt1 = (*m_exp)[i]->GetNumPoints(1);
int pt1 = (*m_exp)[eid]->GetNumPoints(1); int pt2 = (*m_exp)[i]->GetNumPoints(2);
int pt2 = (*m_exp)[eid]->GetNumPoints(2);
int npt0 = (int) pt0*scale; int npt0 = (int) pt0*scale;
int npt1 = (int) pt1*scale; int npt1 = (int) pt1*scale;
int npt2 = (int) pt2*scale; int npt2 = (int) pt2*scale;
LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[eid]->GetPointsType(0)); LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0));
LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[eid]->GetPointsType(1)); LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1));
LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[eid]->GetPointsType(2)); LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2));
// Interpolate points; // Interpolate points;
LibUtilities::Interp3D((*m_exp)[eid]->GetBasis(0)->GetPointsKey(), LibUtilities::Interp3D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(1)->GetPointsKey(), (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(2)->GetPointsKey(), (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
&inarray[cnt], newPointsKey0, &inarray[cnt], newPointsKey0,
newPointsKey1, newPointsKey2, newPointsKey1, newPointsKey2,
&outarray[cnt1]); &outarray[cnt1]);
...@@ -607,26 +606,25 @@ namespace Nektar ...@@ -607,26 +606,25 @@ namespace Nektar
for(int i = 0; i < GetExpSize(); ++i) for(int i = 0; i < GetExpSize(); ++i)
{ {
// get new points key // get new points key
int eid = i; int pt0 = (*m_exp)[i]->GetNumPoints(0);
int pt0 = (*m_exp)[eid]->GetNumPoints(0); int pt1 = (*m_exp)[i]->GetNumPoints(1);
int pt1 = (*m_exp)[eid]->GetNumPoints(1); int pt2 = (*m_exp)[i]->GetNumPoints(2);
int pt2 = (*m_exp)[eid]->GetNumPoints(2);
int npt0 = (int) pt0*scale; int npt0 = (int) pt0*scale;
int npt1 = (int) pt1*scale; int npt1 = (int) pt1*scale;
int npt2 = (int) pt2*scale; int npt2 = (int) pt2*scale;
LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[eid]->GetPointsType(0)); LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0));
LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[eid]->GetPointsType(1)); LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1));
LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[eid]->GetPointsType(2)); LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2));
// Project points; // Project points;
LibUtilities::PhysGalerkinProject3D(newPointsKey0, LibUtilities::PhysGalerkinProject3D(newPointsKey0,
newPointsKey1, newPointsKey1,
newPointsKey2, newPointsKey2,
&inarray[cnt], &inarray[cnt],
(*m_exp)[eid]->GetBasis(0)->GetPointsKey(), (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(1)->GetPointsKey(), (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
(*m_exp)[eid]->GetBasis(2)->GetPointsKey(), (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
&outarray[cnt1]); &outarray[cnt1]);
cnt += npt0*npt1*npt2; cnt += npt0*npt1*npt2;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment