Commit 1670d7df authored by Kilian Lackhove's avatar Kilian Lackhove
Browse files

EquationSystem: Moved pts file import and interpolation logic into separate method

parent c6205f52
......@@ -1000,70 +1000,77 @@ namespace Nektar
}
LibUtilities::PtsFieldSharedPtr outPts;
std::string funcFilename = m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
InterpPts(funcFilename, filename, outPts);
LibUtilities::PtsFieldSharedPtr ptsField;
int fieldInd;
vector<string> fieldNames = outPts->GetFieldNames();
for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
{
if (outPts->GetFieldName(fieldInd) == pFieldName)
{
break;
}
}
ASSERTL0(fieldInd != fieldNames.size(), "field not found");
pArray = outPts->GetPts(fieldInd + outPts->GetDim());
}
void EquationSystem::InterpPts(std::string funcFilename,
std::string filename,
LibUtilities::PtsFieldSharedPtr &outPts)
{
unsigned int nq = m_fields[0]->GetNpoints();
LibUtilities::PtsFieldSharedPtr inPts;
LibUtilities::PtsIO ptsIO(m_session->GetComm());
ptsIO.Import(filename, ptsField);
ptsIO.Import(filename, inPts);
Array<OneD, Array<OneD, NekDouble> > pts(ptsField->GetDim() +
ptsField->GetNFields());
for (int i = 0; i < ptsField->GetDim() + ptsField->GetNFields(); ++i)
Array<OneD, Array<OneD, NekDouble> > pts(inPts->GetDim() + inPts->GetNFields());
for (int i = 0; i < inPts->GetDim() + inPts->GetNFields(); ++i)
{
pts[i] = Array<OneD, NekDouble>(nq);
}
if (ptsField->GetDim() == 1)
if (inPts->GetDim() == 1)
{
m_fields[0]->GetCoords(pts[0]);
}
else if (ptsField->GetDim() == 2)
else if (inPts->GetDim() == 2)
{
m_fields[0]->GetCoords(pts[0], pts[1]);
}
else if (ptsField->GetDim() == 3)
else if (inPts->GetDim() == 3)
{
m_fields[0]->GetCoords(pts[0], pts[1], pts[2]);
}
outPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
ptsField->GetDim(), ptsField->GetFieldNames(), pts);
inPts->GetDim(), inPts->GetFieldNames(), pts);
// check if we already computed this funcKey combination
std::string interpKey =
m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
map<std::string, FieldUtils::Interpolator>::iterator it =
m_interpolators.find(interpKey);
m_interpolators.find(funcFilename);
if (it == m_interpolators.end())
{
m_interpolators[interpKey] =
m_interpolators[funcFilename] =
FieldUtils::Interpolator(Nektar::FieldUtils::eShepard);
if (m_comm->GetRank() == 0)
{
m_interpolators[interpKey].SetProgressCallback(
m_interpolators[funcFilename].SetProgressCallback(
&EquationSystem::PrintProgressbar, this);
}
m_interpolators[interpKey].CalcWeights(ptsField, outPts);
m_interpolators[funcFilename].CalcWeights(inPts, outPts);
if (m_comm->GetRank() == 0)
{
cout << endl;
if (GetSession()->DefinesCmdLineArgument("verbose"))
{
m_interpolators[interpKey].PrintStatistics();
m_interpolators[funcFilename].PrintStatistics();
}
}
}
m_interpolators[interpKey].Interpolate(ptsField, outPts);
int fieldInd;
vector<string> fieldNames = ptsField->GetFieldNames();
for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
{
if (ptsField->GetFieldName(fieldInd) == pFieldName)
{
break;
}
}
ASSERTL0(fieldInd != fieldNames.size(), "field not found");
m_interpolators[funcFilename].Interpolate(inPts, outPts);
pArray = pts[ptsField->GetDim() + fieldInd];
}
......
......@@ -190,6 +190,9 @@ namespace Nektar
const std::string& pFunctionName,
const NekDouble& pTime = 0.0,
const int domain = 0);
SOLVER_UTILS_EXPORT void InterpPts(
std::string funcFilename, std::string filename, Nektar::LibUtilities::PtsFieldSharedPtr& outPts );
// Describe a function.
SOLVER_UTILS_EXPORT std::string DescribeFunction(
......
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