Commit d9397fc8 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Merge remote-tracking branch 'ssherw/fix/Averagefield_varP' into feature/IncNS_Sensor

parents fc48ccc9 3a9cb494
......@@ -20,6 +20,8 @@ v4.4.1
- 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 BLAS CMake dependencies (!763)
- Fix interpolation issue with Lagrange basis functions (!768)
**FieldConvert**:
- Fix issue with FieldConvert when range flag used (!761)
......@@ -29,6 +31,7 @@ v4.4.1
- Rework meshing control so that if possible viewable meshes will be dumped
when some part of the system fails (!756)
- Add manifold meshing option (!756)
- Fix issue with older rea input files (!765)
**FieldConvert:**
- Fix issue with field ordering in the interppointdatatofld module (!754)
......
......@@ -484,6 +484,9 @@ IF( NEKTAR_USE_BLAS_LAPACK )
TARGET_LINK_LIBRARIES(LibUtilities LINK_PUBLIC ${NATIVE_LAPACK} ${NATIVE_BLAS})
ENDIF( NEKTAR_USE_SYSTEM_BLAS_LAPACK )
IF(THIRDPARTY_BUILD_BLAS_LAPACK)
ADD_DEPENDENCIES(LibUtilities lapack-3.7.0)
ENDIF()
ENDIF( NEKTAR_USE_BLAS_LAPACK )
IF( NEKTAR_USE_PETSC )
......
......@@ -1103,9 +1103,7 @@ namespace Nektar
"Not set up for non boundary-interior expansions");
ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
"Assuming that input matrix was square");
ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A,
"Method set up only for modified modal bases curretly");
int i,j;
int id1,id2;
Expansion2DSharedPtr faceExp = m_faceExp[face].lock();
......
......@@ -710,6 +710,26 @@ namespace Nektar
}
break;
}
case LibUtilities::eGLL_Lagrange:
{
LibUtilities::PointsKey
p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey
p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey
p2(nummodes[2], LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey t0(
m_base[0]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey t1(
m_base[1]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey t2(
m_base[2]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::Interp3D(p0, p1, p2, data, t0, t1, t2, coeffs);
}
break;
default:
ASSERTL0(false, "basis is either not set up or not "
"hierarchicial");
......
......@@ -1563,15 +1563,17 @@ namespace Nektar
break;
case LibUtilities::eGLL_Lagrange:
{
// Assume that input is also Gll_Lagrange but no way to check;
LibUtilities::PointsKey
p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey
p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
LibUtilities::Interp2D(p0, p1, data,
m_base[0]->GetPointsKey(),
m_base[1]->GetPointsKey(),
coeffs);
LibUtilities::PointsKey t0(
m_base[0]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey t1(
m_base[1]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
}
break;
case LibUtilities::eGauss_Lagrange:
......@@ -1581,10 +1583,13 @@ namespace Nektar
p0(nummodes[0],LibUtilities::eGaussGaussLegendre);
LibUtilities::PointsKey
p1(nummodes[1],LibUtilities::eGaussGaussLegendre);
LibUtilities::Interp2D(p0, p1, data,
m_base[0]->GetPointsKey(),
m_base[1]->GetPointsKey(),
coeffs);
LibUtilities::PointsKey t0(
m_base[0]->GetNumModes(),
LibUtilities::eGaussGaussLegendre);
LibUtilities::PointsKey t1(
m_base[1]->GetNumModes(),
LibUtilities::eGaussGaussLegendre);
LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
}
break;
default:
......
......@@ -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 = 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);
......@@ -2124,7 +2122,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;
......@@ -2135,20 +2133,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<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
......@@ -2156,19 +2152,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:
......@@ -2177,7 +2173,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;
......@@ -2187,7 +2183,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;
......@@ -2198,29 +2194,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
......@@ -2235,7 +2231,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]);
}
}
......
......@@ -2274,8 +2274,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;
......@@ -2289,30 +2289,29 @@ using namespace boost::assign;
LocalRegions::Expansion3DSharedPtr exp =
(*m_exp)[i]->as<LocalRegions::Expansion3D>();
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<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]);
......@@ -2320,7 +2319,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:
......@@ -2331,7 +2330,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;
......@@ -2343,7 +2342,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;
......@@ -2355,7 +2354,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;
......@@ -2366,34 +2365,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
......@@ -2408,8 +2407,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]);
}
}
......
......@@ -2371,18 +2371,17 @@ namespace Nektar
{
std::vector<unsigned int> nummodes;
vector<LibUtilities::BasisType> 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();
}
}
......
......@@ -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;
......
......@@ -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));