Add support for zlib-encoded input/output, fix to 'official' spacings used in...

Add support for zlib-encoded input/output, fix to 'official' spacings used in boundary condition and curve output

utilities/NekMesh/CMakeLists.txt

@@ -119,6 +119,9 @@ ADD_NEKTAR_TEST        (Gmsh/SquareTri)
 # Nektar tests
 ADD_NEKTAR_TEST        (Nektar/HexLinear)
 ADD_NEKTAR_TEST        (Nektar/Tube45)
+# Nek5000 tests
+ADD_NEKTAR_TEST        (Nek5000/3da)
+ADD_NEKTAR_TEST        (Nek5000/r1854a)
 # StarCCM tests
 ADD_NEKTAR_TEST        (StarTec/CubePer)
 ADD_NEKTAR_TEST_LENGTHY(StarTec/StraightRW)

utilities/NekMesh/InputModules/InputGmsh.cpp

@@ -339,7 +339,7 @@ void InputGmsh::Process()
             }
         }
     }
-    m_mshFile.close();
+    m_mshFile.reset();
 
     // Go through element and remap tags if necessary.
     map<int, map<LibUtilities::ShapeType, int> > compMap;

utilities/NekMesh/InputModules/InputNek.cpp

@@ -896,7 +896,7 @@ void InputNek::Process()
         nSurfaces++;
     }
 
-    m_mshFile.close();
+    m_mshFile.reset();
 
     // -- Process rest of mesh.
     ProcessEdges();

utilities/NekMesh/InputModules/InputNek5000.cpp

@@ -65,13 +65,13 @@ InputNek5000::~InputNek5000()
 }
 
 /**
- * @brief Processes Nektar file format.
+ * @brief Processes Nek5000 file format.
  *
- * Nektar sessions are defined by rea files, and contain sections defining a DNS
- * simulation in a specific order. The converter only reads mesh information,
- * curve information if it exists and boundary information.
- *
- * @param pFilename Filename of Nektar session file to read.
+ * Nek5000 sessions are defined by rea files, and contain sections defining a
+ * DNS simulation in a specific order. The converter only reads mesh
+ * information, curve information if it exists and boundary information. The
+ * format is similar to the rea format supported by #InputNek, but the layout is
+ * sufficiently different that this module is separate.
  */
 void InputNek5000::Process()
 {
@@ -192,7 +192,7 @@ void InputNek5000::Process()
         }
         else
         {
-            // - hex: 4 lines with x/y/z-coords for base 4 nodes, then 4 more
+            // - hex: 3 lines with x/y/z-coords for base 4 nodes, then 3 more
             //   for upper 4 nodes
             elType = LibUtilities::eHexahedron;
             nNodes = 8;
@@ -209,9 +209,9 @@ void InputNek5000::Process()
             }
         }
 
-        // Nektar meshes do not contain a unique list of nodes, so this
-        // block constructs a unique set so that elements can be created
-        // with unique nodes.
+        // Nek5000 meshes do not contain a unique list of nodes, so this block
+        // constructs a unique set so that elements can be created with unique
+        // nodes.
         vector<NodeSharedPtr> nodeList(nNodes);
         for (k = 0; k < nNodes; ++k)
         {
@@ -276,14 +276,45 @@ void InputNek5000::Process()
         for (i = 0; i < nCurves; ++i)
         {
             getline(m_mshFile, line);
-            s.clear();
-            s.str(line);
 
             int elmt, side;
             NekDouble curveData[5];
             char curveType;
 
-            s >> side >> elmt;
+            if (nElements < 1000)
+            {
+                // side in first 3 characters, elmt in next 3
+                s.str(line.substr(0, 3));
+                s >> side;
+                s.clear();
+                s.str(line.substr(3, 3));
+                s >> elmt;
+                line = line.substr(6);
+            }
+            else if (nElements < 1000000)
+            {
+                // side in first 2 characters, elmt in next 6
+                s.str(line.substr(0, 2));
+                s >> side;
+                s.clear();
+                s.str(line.substr(2, 6));
+                s >> elmt;
+                line = line.substr(8);
+            }
+            else
+            {
+                // side in first 2 characters, elmt in next 12
+                s.str(line.substr(0, 2));
+                s >> side;
+                s.clear();
+                s.str(line.substr(2, 12));
+                s >> elmt;
+                line = line.substr(14);
+            }
+
+            s.clear();
+            s.str(line);
+
             for (j = 0; j < 5; ++j)
             {
                 s >> curveData[j];
@@ -352,7 +383,8 @@ void InputNek5000::Process()
                     link      = centroid - midpoint;
                     sign      = link.dot(unitNormal);
                     sign      = convexity * sign / fabs(sign);
-                    centre    = midpoint + unitNormal * (sign * cos(semiangle) * radius);
+                    centre    = midpoint + unitNormal * (sign * cos(semiangle) *
+                                                         radius);
 
                     NekDouble theta1, theta2, dtheta, phi;
                     theta1 = atan2 (P1.m_y - centre.m_y, P1.m_x - centre.m_x);
@@ -404,14 +436,15 @@ void InputNek5000::Process()
     }
 
     int nSurfaces = 0;
-    boost::unordered_set<pair<int, int> > periodicIn, periodicOut;
+    boost::unordered_set<pair<int, int> > periodicIn;
     int periodicInId = -1, periodicOutId = -1;
 
+    // Boundary conditions: should be precisely nElements * nFaces lines to
+    // read.
+    int lineCnt = 0;
     while (m_mshFile.good())
     {
         getline(m_mshFile, line);
-        s.clear();
-        s.str(line);
 
         // Found a new section. We don't support anything in the rea file beyond
         // this point so we'll just quit.
@@ -425,7 +458,55 @@ void InputNek5000::Process()
         int elmt, side;
         NekDouble data[5];
 
-        s >> bcType >> elmt >> side;
+        // type in chars 0-3
+        s.clear();
+        s.str(line.substr(0, 4));
+        s >> bcType;
+
+        if (nElements < 1000)
+        {
+            // elmt in chars 4-6, side in next 3
+            s.clear();
+            s.str(line.substr(4, 3));
+            s >> elmt;
+            s.clear();
+            s.str(line.substr(7, 3));
+            s >> side;
+            line = line.substr(10);
+        }
+        else if (nElements < 100000)
+        {
+            // elmt in chars 4-8, side in next 1
+            s.clear();
+            s.str(line.substr(4, 5));
+            s >> elmt;
+            s.clear();
+            s.str(line.substr(9, 1));
+            s >> side;
+            line = line.substr(10);
+        }
+        else if (nElements < 1000000)
+        {
+            // elmt in chars 4-9, no side
+            s.clear();
+            s.str(line.substr(4, 6));
+            s >> elmt;
+            side = lineCnt % (2 * m_mesh->m_expDim);
+            line = line.substr(9);
+        }
+        else
+        {
+            // elmt in chars 4-15, no side
+            s.clear();
+            s.str(line.substr(4, 12));
+            s >> elmt;
+            side = lineCnt % (2 * m_mesh->m_expDim);
+            line = line.substr(15);
+        }
+
+        s.clear();
+        s.str(line);
+
         for (i = 0; i < 5; ++i)
         {
             s >> data[i];
@@ -435,6 +516,9 @@ void InputNek5000::Process()
         --elmt;
         --side;
 
+        // Increment lines read
+        lineCnt++;
+
         ElementSharedPtr el = m_mesh->m_element[m_mesh->m_spaceDim][elmt];
 
         std::string fields[] = { "u", "v", "w", "p" };
@@ -611,7 +695,13 @@ void InputNek5000::Process()
         nSurfaces++;
     }
 
-    m_mshFile.close();
+    if (lineCnt != nElements * (m_mesh->m_expDim * 2))
+    {
+        cerr << "Warning: boundary conditions may not have been correctly read "
+             << "from Nek5000 input file." << endl;
+    }
+
+    m_mshFile.reset();
 
     // -- Process rest of mesh.
     ProcessEdges();

utilities/NekMesh/InputModules/InputPly.cpp

@@ -70,7 +70,7 @@ void InputPly::Process()
 
     ReadPly(m_mshFile);
 
-    m_mshFile.close();
+    m_mshFile.reset();
 
     ProcessVertices();
     ProcessEdges();
@@ -79,7 +79,7 @@ void InputPly::Process()
     ProcessComposites();
 }
 
-void InputPly::ReadPly(std::ifstream &mshFile, NekDouble scale)
+void InputPly::ReadPly(io::filtering_istream &mshFile, NekDouble scale)
 {
     m_mesh->m_expDim = 0;
     string line;

utilities/NekMesh/InputModules/InputPly.h

@@ -60,7 +60,7 @@ public:
     /// Populate and validate required data structures.
     virtual void Process();
 
-    void ReadPly(std::ifstream &mshFile, NekDouble scale = 1.0);
+    void ReadPly(io::filtering_istream &mshFile, NekDouble scale = 1.0);
 
 private:
 };

utilities/NekMesh/InputModules/InputSem.cpp

@@ -649,7 +649,7 @@ void InputSem::Process()
     }
 
     PrintSummary();
-    m_mshFile.close();
+    m_mshFile.reset();
 
     // Process rest of mesh.
     ProcessVertices();

utilities/NekMesh/InputModules/InputStarTec.cpp

@@ -109,7 +109,7 @@ void InputTec::Process()
     }
 
     PrintSummary();
-    m_mshFile.close();
+    m_mshFile.reset();
 
     ProcessEdges();
     ProcessFaces();
@@ -125,15 +125,12 @@ void InputTec::ReadZone(int &nComposite)
     int start, end;
     stringstream s;
     NekDouble value;
-    streampos pos;
     static int zcnt = 1;
 
     // Read Zone Header
     nnodes = nfaces = nelements = 0;
     while (!m_mshFile.eof())
     {
-        pos = m_mshFile.tellg();
-
         getline(m_mshFile, line);
 
         boost::to_upper(line);
@@ -141,7 +138,6 @@ void InputTec::ReadZone(int &nComposite)
         // cehck to see if readable data.
         if (sscanf(line.c_str(), "%lf", &value) == 1)
         {
-            m_mshFile.seekg(pos);
             break;
         }
 
@@ -200,31 +196,34 @@ void InputTec::ReadZone(int &nComposite)
 
     cout << "Setting up zone " << zcnt++;
 
-    vector<NekDouble> x, y, z;
+    int nodeCount = 3 * nnodes;
+    vector<NekDouble> nodeLocs;
 
-    // Read in Nodes
-    for (i = 0; i < nnodes; ++i)
+    while (nodeCount > 0 && !m_mshFile.eof())
     {
-        m_mshFile >> value;
-        x.push_back(value);
-    }
-
-    for (i = 0; i < nnodes; ++i)
-    {
-        m_mshFile >> value;
-        y.push_back(value);
+        s.clear();
+        s.str(line);
+        while (s >> value)
+        {
+            nodeLocs.push_back(value);
+            nodeCount--;
+        }
+        if (nodeCount > 0)
+        {
+            getline(m_mshFile, line);
+        }
     }
 
-    for (i = 0; i < nnodes; ++i)
-    {
-        m_mshFile >> value;
-        z.push_back(value);
-    }
+    ASSERTL0(nodeLocs.size() == 3*nnodes, "Unable to read correct number of "
+             "nodes from Tecplot file");
 
     std::vector<NodeSharedPtr> Nodes;
     for (i = 0; i < nnodes; ++i)
     {
-        Nodes.push_back(boost::shared_ptr<Node>(new Node(i, x[i], y[i], z[i])));
+        Nodes.push_back(
+            boost::shared_ptr<Node>(
+                new Node(i, nodeLocs[i], nodeLocs[i+nnodes],
+                         nodeLocs[i+2*nnodes])));
     }
 
     // Read Node count per face

utilities/NekMesh/InputModules/InputSwan.cpp

@@ -84,7 +84,7 @@ void InputSwan::Process()
     if (tmp[0] != tmp[5] || tmp[0] != 4 * sizeof(int))
     {
         cout << "Header data broken" << endl;
-        m_mshFile.close();
+        m_mshFile.reset();
         return;
     }
 
@@ -109,7 +109,7 @@ void InputSwan::Process()
         tets[0] != ND * NB_Tet * sizeof(int))
     {
         cout << "ERROR [InputSwan]: Tetrahedron data broken." << endl;
-        m_mshFile.close();
+        m_mshFile.reset();
         return;
     }
 
@@ -126,7 +126,7 @@ void InputSwan::Process()
     if (tmp[0] != tmp[1] || tmp[0] != 3 * NB_Points * sizeof(double))
     {
         cout << "ERROR [InputSwan]: Point data broken." << endl;
-        m_mshFile.close();
+        m_mshFile.reset();
         return;
     }
 
@@ -180,7 +180,7 @@ void InputSwan::Process()
     if (tmp[0] != tmp[1])
     {
         cout << "ERROR [InputSwan]: Surface data broken." << endl;
-        m_mshFile.close();
+        m_mshFile.reset();
         return;
     }
 
@@ -206,7 +206,7 @@ void InputSwan::Process()
         m_mesh->m_element[2].push_back(E);
     }
 
-    m_mshFile.close();
+    m_mshFile.reset();
 
     // Process the rest of the mesh.
     ProcessVertices();

utilities/NekMesh/Module.cpp

@@ -33,9 +33,12 @@
 //
 ////////////////////////////////////////////////////////////////////////////////
 
+#include <boost/iostreams/filter/gzip.hpp>
+
 #include "Module.h"
 
 using namespace std;
+namespace io = boost::iostreams;
 
 namespace Nektar
 {
@@ -78,7 +81,19 @@ OutputModule::OutputModule(MeshSharedPtr m) : Module(m)
 void InputModule::OpenStream()
 {
     string fname = m_config["infile"].as<string>();
-    m_mshFile.open(fname.c_str());
+
+    if (fname.size() > 3 && fname.substr(fname.size() - 3, 3) == ".gz")
+    {
+        m_mshFileStream.open(fname.c_str(), ios_base::in | ios_base::binary);
+        m_mshFile.push(io::gzip_decompressor());
+        m_mshFile.push(m_mshFileStream);
+    }
+    else
+    {
+        m_mshFileStream.open(fname.c_str());
+        m_mshFile.push(m_mshFileStream);
+    }
+
     if (!m_mshFile.good())
     {
         cerr << "Error opening file: " << fname << endl;
@@ -92,7 +107,19 @@ void InputModule::OpenStream()
 void OutputModule::OpenStream()
 {
     string fname = m_config["outfile"].as<string>();
-    m_mshFile.open(fname.c_str());
+
+    if (fname.size() > 3 && fname.substr(fname.size() - 3, 3) == ".gz")
+    {
+        m_mshFileStream.open(fname.c_str(), ios_base::out | ios_base::binary);
+        m_mshFile.push(io::gzip_compressor());
+        m_mshFile.push(m_mshFileStream);
+    }
+    else
+    {
+        m_mshFileStream.open(fname.c_str());
+        m_mshFile.push(m_mshFileStream);
+    }
+
     if (!m_mshFile.good())
     {
         cerr << "Error opening file: " << fname << endl;

utilities/NekMesh/Module.h

@@ -36,6 +36,9 @@
 #ifndef UTILITIES_NEKMESH_MODULE
 #define UTILITIES_NEKMESH_MODULE
 
+#include <boost/iostreams/device/file_descriptor.hpp>
+#include <boost/iostreams/filtering_stream.hpp>
+
 #include <map>
 #include <vector>
 #include <sstream>
@@ -47,6 +50,8 @@
 #include <LibUtilities/BasicUtils/NekFactory.hpp>
 #include <NekMeshUtils/MeshElements/Mesh.h>
 
+namespace io = boost::iostreams;
+
 namespace Nektar
 {
     namespace Utilities
@@ -209,9 +214,11 @@ namespace Nektar
 
         protected:
             /// Print summary of elements.
-            void         PrintSummary();
+            void PrintSummary();
             /// Input stream
-            std::ifstream m_mshFile;
+            io::filtering_istream m_mshFile;
+            /// Input stream
+            std::ifstream m_mshFileStream;
         };
 
         /**
@@ -243,7 +250,9 @@ namespace Nektar
 
         protected:
             /// Output stream
-            std::ofstream m_mshFile;
+            io::filtering_ostream m_mshFile;
+            /// Input stream
+            std::ofstream m_mshFileStream;
         };
 
         typedef std::pair<ModuleType,std::string> ModuleKey;

utilities/NekMesh/NekMesh.cpp

@@ -200,6 +200,14 @@ int main(int argc, char* argv[])
             {
                 int    dot    = tmp1[0].find_last_of('.') + 1;
                 string ext    = tmp1[0].substr(dot, tmp1[0].length() - dot);
+
+                if (ext == "gz")
+                {
+                    string tmp = tmp1[0].substr(0, tmp1[0].length() - 3);
+                    dot = tmp.find_last_of('.') + 1;
+                    ext = tmp.substr(dot, tmp.length() - dot);
+                }
+
                 module.second = ext;
                 tmp1.push_back(string(i == 0 ? "infile=" : "outfile=")+tmp1[0]);
             }

utilities/NekMesh/ProcessModules/ProcessSpherigon.cpp

@@ -474,11 +474,15 @@ void ProcessSpherigon::Process()
                  << " with scaling of " << scale << endl;
         }
 
-        ifstream inply;
+        ifstream inplyTmp;
+        io::filtering_istream inply;
         InputPlySharedPtr plyfile;
 
-        inply.open(normalfile.c_str());
-        ASSERTL0(inply, string("Could not open input ply file: ") + normalfile);
+        inplyTmp.open(normalfile.c_str());
+        ASSERTL0(inplyTmp,
+                 string("Could not open input ply file: ") + normalfile);
+
+        inply.push(inplyTmp);
 
         int j;
         MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());