fix interppoints module

library/FieldUtils/ProcessModules/ProcessInterpPoints.cpp

@@ -45,6 +45,8 @@ using namespace std;
 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
 #include <LibUtilities/BasicUtils/Progressbar.hpp>
 #include <LibUtilities/BasicUtils/SharedArray.hpp>
+#include <LibUtilities/BasicUtils/PtsIO.h>
+#include <LibUtilities/BasicUtils/CsvIO.h>
 #include <boost/lexical_cast.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 
@@ -116,6 +118,12 @@ void ProcessInterpPoints::Process(po::variables_map &vm)
         MemoryManager<SpatialDomains::DomainRange>::AllocateSharedPtr();
     int coordim = m_f->m_fieldPts->GetDim();
     int npts    = m_f->m_fieldPts->GetNpoints();
+    std::vector<std::string> fieldNames = m_f->m_fieldPts->GetFieldNames();
+    for (auto it = fieldNames.begin(); it != fieldNames.end(); ++it)
+    {
+        m_f->m_fieldPts->RemoveField(*it);
+    }
+
     Array<OneD, Array<OneD, NekDouble> > pts;
     m_f->m_fieldPts->GetPts(pts);
     rng->m_checkShape = false;
@@ -216,14 +224,30 @@ void ProcessInterpPoints::CreateFieldPts(po::variables_map &vm)
     int rank   = m_f->m_comm->GetRank();
     int nprocs = m_f->m_comm->GetSize();
     // Check for command line point specification
+
     if (m_config["topts"].as<string>().compare("NotSet") != 0)
     {
         string inFile = m_config["topts"].as<string>();
 
-        LibUtilities::PtsIOSharedPtr ptsIO =
-            MemoryManager<LibUtilities::PtsIO>::AllocateSharedPtr(m_f->m_comm);
+        if (boost::filesystem::path(inFile).extension() == ".pts")
+        {
+            LibUtilities::PtsIOSharedPtr ptsIO =
+                MemoryManager<LibUtilities::PtsIO>::AllocateSharedPtr(m_f->m_comm);
+
+            ptsIO->Import(inFile, m_f->m_fieldPts);
+        }
+        else if (boost::filesystem::path(inFile).extension() == ".csv")
+        {
+            LibUtilities::CsvIOSharedPtr csvIO =
+                MemoryManager<LibUtilities::CsvIO>::AllocateSharedPtr(m_f->m_comm);
+
+            csvIO->Import(inFile, m_f->m_fieldPts);
+        }
+        else
+        {
+            ASSERTL0(false, "no input file found");
+        }
 
-        ptsIO->Import(inFile, m_f->m_fieldPts);
     }
     else if (m_config["line"].as<string>().compare("NotSet") != 0)
     {
@@ -432,7 +456,7 @@ void ProcessInterpPoints::InterpolateFieldToPts(
     NekDouble clamp_up,
     NekDouble def_value)
 {
-    ASSERTL0(pts->GetNFields() >= field0.size(), "ptField has too few fields");
+    ASSERTL0(pts->GetNFields() == field0.size(), "ptField has too few fields");
 
     int nfields = field0.size();
 

library/LibUtilities/BasicUtils/PtsField.cpp

@@ -138,6 +138,26 @@ void PtsField::AddField(const Array< OneD, NekDouble > &pts,
     m_fieldNames.push_back(fieldName);
 }
 
+
+void PtsField::RemoveField(const string fieldName)
+{
+    int nTotvars = m_pts.num_elements();
+
+    // redirect existing pts
+    Array<OneD, Array<OneD, NekDouble> > newpts(nTotvars - 1);
+    for (int i = 0, j = 0; i < nTotvars; ++i)
+    {
+        if (i < GetDim() || m_fieldNames[i - GetDim()] != fieldName)
+        {
+            newpts[j++] = m_pts[i];
+        }
+    }
+
+    m_pts = newpts;
+
+    m_fieldNames.erase(remove(m_fieldNames.begin(), m_fieldNames.end(), fieldName), m_fieldNames.end());
+}
+
 void PtsField::AddPoints(const Array<OneD, const Array<OneD, NekDouble> > &pts)
 {
     ASSERTL1(pts.num_elements() == m_pts.num_elements(),

library/LibUtilities/BasicUtils/PtsField.h

@@ -115,6 +115,8 @@ public:
     LIB_UTILITIES_EXPORT void AddField(const Array<OneD, NekDouble> &pts,
                                        const std::string fieldName);
 
+    LIB_UTILITIES_EXPORT void RemoveField(const std::string fieldName);
+
     LIB_UTILITIES_EXPORT void AddPoints(const Array< OneD, const Array< OneD, NekDouble > > &pts);
 
     LIB_UTILITIES_EXPORT int GetNpoints() const;