diff --git a/library/MultiRegions/DisContField2D.cpp b/library/MultiRegions/DisContField2D.cpp index 885a6212dd124b635d61715642b86c1179ff9e84..1a8c48511843b39439725a84dd37511d92be1922 100644 --- a/library/MultiRegions/DisContField2D.cpp +++ b/library/MultiRegions/DisContField2D.cpp @@ -1814,7 +1814,7 @@ namespace Nektar StdRegions::Orientation edgedir; - int eid,cnt; + int cnt; int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs(); Array 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 = 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(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); @@ -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 loc_lambda(LocBndCoeffs), edge_lambda; Array tmp_coeffs; @@ -2134,20 +2132,18 @@ namespace Nektar // Calculate Q using standard DG formulation. for(i = cnt = 0; i < GetExpSize(); ++i) { - eid = 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(nq_elmt); qrhs1 = Array(nq_elmt); force = Array(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 > 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(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((*m_exp)[eid]->GetGeom()); + SpatialDomains::QuadGeomSharedPtr qGeom = boost::dynamic_pointer_cast((*m_exp)[i]->GetGeom()); ppExp = MemoryManager::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((*m_exp)[eid]->GetGeom()); + SpatialDomains::TriGeomSharedPtr tGeom = boost::dynamic_pointer_cast((*m_exp)[i]->GetGeom()); ppExp = MemoryManager::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 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]); } } diff --git a/library/MultiRegions/ExpList.cpp b/library/MultiRegions/ExpList.cpp index a6f068757906c8380125746d66eef3d43c2d0013..b2c0d1643fd2425a3cd44040db76e665b95162cd 100644 --- a/library/MultiRegions/ExpList.cpp +++ b/library/MultiRegions/ExpList.cpp @@ -2359,18 +2359,17 @@ namespace Nektar { std::vector nummodes; vector basisTypes; - int eid = i; - for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j) + for(int j= 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j) { - nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j)); - basisTypes.push_back(fromExpList->GetExp(eid)->GetBasisType(j)); + nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j)); + basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j)); } - (*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0, - &toCoeffs[m_coeff_offset[eid]], + (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0, + &toCoeffs[m_coeff_offset[i]], basisTypes); - offset += fromExpList->GetExp(eid)->GetNcoeffs(); + offset += fromExpList->GetExp(i)->GetNcoeffs(); } } diff --git a/library/MultiRegions/ExpList2D.cpp b/library/MultiRegions/ExpList2D.cpp index 3d864864d5f9fc2d717421e3c835e4c447419bbe..80630eaf88653fba55fad6cb3b31aaf8265293c7 100644 --- a/library/MultiRegions/ExpList2D.cpp +++ b/library/MultiRegions/ExpList2D.cpp @@ -1227,20 +1227,19 @@ namespace Nektar for(int i = 0; i < GetExpSize(); ++i) { // get new points key - int eid = i; - int pt0 = (*m_exp)[eid]->GetNumPoints(0); - int pt1 = (*m_exp)[eid]->GetNumPoints(1); + int pt0 = (*m_exp)[i]->GetNumPoints(0); + int pt1 = (*m_exp)[i]->GetNumPoints(1); int npt0 = (int) pt0*scale; int npt1 = (int) pt1*scale; LibUtilities::PointsKey newPointsKey0(npt0, - (*m_exp)[eid]->GetPointsType(0)); + (*m_exp)[i]->GetPointsType(0)); LibUtilities::PointsKey newPointsKey1(npt1, - (*m_exp)[eid]->GetPointsType(1)); + (*m_exp)[i]->GetPointsType(1)); // Interpolate points; - LibUtilities::Interp2D((*m_exp)[eid]->GetBasis(0)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(1)->GetPointsKey(), + LibUtilities::Interp2D((*m_exp)[i]->GetBasis(0)->GetPointsKey(), + (*m_exp)[i]->GetBasis(1)->GetPointsKey(), &inarray[cnt],newPointsKey0, newPointsKey1,&outarray[cnt1]); @@ -1260,24 +1259,23 @@ namespace Nektar for(int i = 0; i < GetExpSize(); ++i) { // get new points key - int eid = i; - int pt0 = (*m_exp)[eid]->GetNumPoints(0); - int pt1 = (*m_exp)[eid]->GetNumPoints(1); + int pt0 = (*m_exp)[i]->GetNumPoints(0); + int pt1 = (*m_exp)[i]->GetNumPoints(1); int npt0 = (int) pt0*scale; int npt1 = (int) pt1*scale; LibUtilities::PointsKey newPointsKey0(npt0, - (*m_exp)[eid]->GetPointsType(0)); + (*m_exp)[i]->GetPointsType(0)); LibUtilities::PointsKey newPointsKey1(npt1, - (*m_exp)[eid]->GetPointsType(1)); + (*m_exp)[i]->GetPointsType(1)); // Project points; LibUtilities::PhysGalerkinProject2D( newPointsKey0, newPointsKey1, &inarray[cnt], - (*m_exp)[eid]->GetBasis(0)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(1)->GetPointsKey(), + (*m_exp)[i]->GetBasis(0)->GetPointsKey(), + (*m_exp)[i]->GetBasis(1)->GetPointsKey(), &outarray[cnt1]); cnt += npt0*npt1; diff --git a/library/MultiRegions/ExpList3D.cpp b/library/MultiRegions/ExpList3D.cpp index 4d8e6ef858c904332386ce84fab11cfc5c73869c..b1f8f380db25e87a24159b0831f9b1fc8bc2a5c9 100644 --- a/library/MultiRegions/ExpList3D.cpp +++ b/library/MultiRegions/ExpList3D.cpp @@ -572,22 +572,21 @@ namespace Nektar for(int i = 0; i < GetExpSize(); ++i) { // get new points key - int eid = i; - int pt0 = (*m_exp)[eid]->GetNumPoints(0); - int pt1 = (*m_exp)[eid]->GetNumPoints(1); - int pt2 = (*m_exp)[eid]->GetNumPoints(2); + int pt0 = (*m_exp)[i]->GetNumPoints(0); + int pt1 = (*m_exp)[i]->GetNumPoints(1); + int pt2 = (*m_exp)[i]->GetNumPoints(2); int npt0 = (int) pt0*scale; int npt1 = (int) pt1*scale; int npt2 = (int) pt2*scale; - LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[eid]->GetPointsType(0)); - LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[eid]->GetPointsType(1)); - LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[eid]->GetPointsType(2)); + LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0)); + LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1)); + LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2)); // Interpolate points; - LibUtilities::Interp3D((*m_exp)[eid]->GetBasis(0)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(1)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(2)->GetPointsKey(), + LibUtilities::Interp3D((*m_exp)[i]->GetBasis(0)->GetPointsKey(), + (*m_exp)[i]->GetBasis(1)->GetPointsKey(), + (*m_exp)[i]->GetBasis(2)->GetPointsKey(), &inarray[cnt], newPointsKey0, newPointsKey1, newPointsKey2, &outarray[cnt1]); @@ -607,26 +606,25 @@ namespace Nektar for(int i = 0; i < GetExpSize(); ++i) { // get new points key - int eid = i; - int pt0 = (*m_exp)[eid]->GetNumPoints(0); - int pt1 = (*m_exp)[eid]->GetNumPoints(1); - int pt2 = (*m_exp)[eid]->GetNumPoints(2); + int pt0 = (*m_exp)[i]->GetNumPoints(0); + int pt1 = (*m_exp)[i]->GetNumPoints(1); + int pt2 = (*m_exp)[i]->GetNumPoints(2); int npt0 = (int) pt0*scale; int npt1 = (int) pt1*scale; int npt2 = (int) pt2*scale; - LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[eid]->GetPointsType(0)); - LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[eid]->GetPointsType(1)); - LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[eid]->GetPointsType(2)); + LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0)); + LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1)); + LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2)); // Project points; LibUtilities::PhysGalerkinProject3D(newPointsKey0, newPointsKey1, newPointsKey2, &inarray[cnt], - (*m_exp)[eid]->GetBasis(0)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(1)->GetPointsKey(), - (*m_exp)[eid]->GetBasis(2)->GetPointsKey(), + (*m_exp)[i]->GetBasis(0)->GetPointsKey(), + (*m_exp)[i]->GetBasis(1)->GetPointsKey(), + (*m_exp)[i]->GetBasis(2)->GetPointsKey(), &outarray[cnt1]); cnt += npt0*npt1*npt2;