diff --git a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index 6afe83d8cf78a9f347e7b82ee863b828d0f70b49..20bf6e11efffd3cfa29b83aa835f21567a000277 100644
--- a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -32,7 +32,7 @@
 #include "imstkInternalForceModel.h"
 #include "imstkForceModelConfig.h"
 #include "imstkNonlinearSystem.h"
-#include "imstkVegaMeshReader.h"
+#include "imstkVegaMeshIO.h"
 #include "imstkNewtonSolver.h"
 
 //force models
@@ -40,7 +40,7 @@
 #include "imstkLinearFEMForceModel.h"
 #include "imstkCorotationalFEMForceModel.h"
 #include "imstkIsotropicHyperelasticFEMForceModel.h"
-#include "imstkVegaMeshReader.h"
+#include "imstkVegaMeshIO.h"
 
 //vega
 #include "sparseMatrix.h"
diff --git a/Base/Geometry/Reader/imstkMSHMeshIO.cpp b/Base/Geometry/Reader/imstkMSHMeshIO.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb0e4c65b29e7b143a603139b34eaf0f9092a824
--- /dev/null
+++ b/Base/Geometry/Reader/imstkMSHMeshIO.cpp
@@ -0,0 +1,457 @@
+/*=========================================================================
+
+    Library: iMSTK
+
+    Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+    & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0.txt
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+
+    =========================================================================*/
+
+#include <fstream>
+#include <iostream>
+
+#include "imstkMSHMeshIO.h"
+#include "imstkTetrahedralMesh.h"
+
+#include "tetMesh.h"
+
+#include "g3log/g3log.hpp"
+
+namespace imstk
+{
+
+std::shared_ptr<imstk::VolumetricMesh>
+MSHMeshIO::read(const std::string & filePath, const MeshFileType meshType)
+{
+    if (meshType != MeshFileType::MSH)
+    {
+        LOG(WARNING) <<"MSHMeshIO::read error: file type not supported";
+        return nullptr;
+    }
+
+    // based on the format provided on
+    // http://www.manpagez.com/info/gmsh/gmsh-2.2.6/gmsh_63.php
+
+    size_t Nnodes; // number of nodes
+    std::vector<size_t> Nodes_tag;                                                 // number assigned to each node (node number)
+    std::vector<double> Nodes_xCoor, Nodes_yCoor, Nodes_zCoor;                     // nodes coordinates
+    size_t N1Delems = 0;                                                           // total number of lines
+    size_t N2Delems = 0;                                                           // total number of surface elements (considering only triangle for now), cannot handle mix surface element types (ex. tri+quad, and so on).
+    size_t N3Delems = 0;                                                           // total number of volumetric elements  (considering only tets for now)
+    std::vector<size_t> N1Delems_tag;                                              // number assigned to each line element
+    std::vector<size_t> N2Delems_tag;                                              // number assigned to each surface (triangle) element
+    std::vector<size_t> N3Delems_tag;                                              // number assigned to each volumetric (tet) element
+    std::vector<std::array<size_t,2>> N1Delems_conn;                               // line element connectivity
+    std::vector<std::array<size_t,3>> N2Delems_conn;                               // triangle element connectivity
+    std::vector<std::array<size_t,4>> N3Delems_conn;                               // tet element connectivity
+    size_t OtherElems;                                                             // all element types other than line, triangle, tets.
+    std::vector<int> tag_positions;
+
+    std::ifstream msh_stream(filePath);
+    if (!msh_stream.is_open())
+    {
+        LOG(WARNING) << "Failed to open the input .msh file";
+        return nullptr;
+    }
+
+    std::string msh_line;
+    std::stringstream msh_ss;
+    std::string key_word;
+
+    // look for "$MeshFormat"
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+            if (key_word == "$MeshFormat")
+            {
+                break;
+            }
+        }
+    }
+
+    if (msh_stream.eof())
+    {
+        LOG(INFO) << "Warning:  version number, file-type, data-size not found in the msh file.";
+    }
+
+    msh_stream.clear();
+    msh_stream.seekg(0, std::ios::beg);
+
+    // look for "$NodeData"
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+
+            if (!key_word.compare("$NOD") || !key_word.compare("$Nodes"))
+            {
+                break;
+            }
+        }
+    }
+
+    if (msh_stream.eof())
+    {
+        LOG(WARNING) << "Error: Nodes not defined." ;
+        return nullptr;
+    }
+    else
+    {
+        std::string nnodes;
+        while (getline(msh_stream, msh_line))
+        {
+            msh_ss.str(std::string());
+            msh_ss.clear();
+            if (!msh_line.empty())
+            {
+                msh_ss << msh_line;
+                msh_ss >> nnodes;
+                Nnodes = stoi(nnodes);
+                break;
+            }
+        }
+    }
+
+    LOG(INFO) << "The MSH mesh comprises of: \n" << '\t' << "Number of NODES: " << Nnodes ;
+
+    // get coordinates ( geometry )
+
+    Nodes_tag.resize(Nnodes);
+    Nodes_xCoor.resize(Nnodes);
+    Nodes_yCoor.resize(Nnodes);
+    Nodes_zCoor.resize(Nnodes);
+
+    std::string node_tag;
+    std::string node_xC;
+    std::string node_yC;
+    std::string node_zC;
+    size_t nodes_count = 0;
+
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+            if (!key_word.compare("$ENDNOD") || !key_word.compare("$EndNodes"))
+            {
+                break;
+            }
+            else
+            {
+                msh_ss.str(std::string());
+                msh_ss.clear();
+                msh_ss << msh_line;
+                msh_ss >> node_tag;
+                Nodes_tag[nodes_count] = stoul(node_tag);
+                // x coordinate
+                msh_ss >> node_xC;
+                Nodes_xCoor[nodes_count] = stod(node_xC);
+                // y coordinate
+                msh_ss >> node_yC;
+                Nodes_yCoor[nodes_count] = stod(node_yC);
+                // z coordinate
+                msh_ss >> node_zC;
+                Nodes_zCoor[nodes_count] = stod(node_zC);
+                ++nodes_count;
+            }
+        }
+    }
+
+    if (nodes_count != Nnodes)
+    {
+        LOG(WARNING) << "Error in reading the nodes from the input MSH file.";
+        return nullptr;
+    }
+
+    // Get the elements ( topology )
+    msh_stream.clear();
+    msh_stream.seekg(0, std::ios::beg);
+
+    // look for "$ELM"
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+            if (!key_word.compare("$ELM") || !key_word.compare("$Elements"))
+            {
+                break;
+            }
+        }
+    }
+
+    if (msh_stream.eof())
+    {
+        LOG(WARNING) << "Error: Elements not defined.";
+        return nullptr;
+    }
+    else
+    {
+        // get the total number of elements of each topology type
+        size_t lineEl_count = 0;
+        size_t surfEl_count = 0;
+        size_t volEl_count = 0;
+        size_t otherEl_count = 0;
+        size_t elem_type;
+        std::string temp;
+
+        while (getline(msh_stream, msh_line))
+        {
+            msh_ss.str(std::string());
+            msh_ss.clear();
+            if (!msh_line.empty())
+            {
+                break;
+            }
+        }
+
+        while (getline(msh_stream, msh_line))
+        {
+            msh_ss.str(std::string());
+            msh_ss.clear();
+            if (!msh_line.empty())
+            {
+                msh_ss << msh_line;
+                msh_ss >> key_word;
+                if (!key_word.compare("$ENDELM") || !key_word.compare("$EndElements"))
+                {
+                    break;
+                }
+                else
+                {
+                    msh_ss.str(std::string());
+                    msh_ss.clear();
+                    msh_ss << msh_line;
+                    msh_ss >> temp;
+                    msh_ss >> temp;
+                    elem_type = stoul(temp);
+                    if (elem_type == 1)
+                    {
+                        ++lineEl_count;
+                    }
+                    else if (elem_type == 2)
+                    {
+                        ++surfEl_count;
+                    }
+                    else if (elem_type == 3)
+                    {
+                        ++otherEl_count;
+                    }
+                    else if (elem_type == 4)
+                    {
+                        ++volEl_count;
+                    }
+                    else if (elem_type >= 5 && elem_type <= 32)
+                    {
+                        ++otherEl_count;
+                    }
+                    else
+                    {
+                        LOG(WARNING) << "Specified wrong element types.";
+                        return nullptr;
+                    }
+                }
+            }
+        }
+        N1Delems = lineEl_count;
+        N2Delems = surfEl_count;
+        N3Delems = volEl_count;
+        OtherElems = otherEl_count;
+    }
+
+    // set the stream back to the elem field
+    msh_stream.clear();
+    msh_stream.seekg(0, std::ios::beg);
+
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+            if (!key_word.compare("$ELM") || !key_word.compare("$Elements"))
+            {
+                break;
+            }
+        }
+    }
+
+    getline(msh_stream, msh_line);
+    std::string temp;
+    msh_ss.str(std::string());
+    msh_ss.clear();
+    msh_ss << msh_line;
+    msh_ss >> temp;
+
+    if (!(stoul(temp) == (N1Delems + N2Delems + N3Delems + OtherElems)))
+    {
+        LOG(WARNING) << "Error reading the element field in the msh file .. exiting";
+        return nullptr;
+    }
+    if (N3Delems == 0)
+    {
+        LOG(WARNING) << "No volumetric ( tetrahedral element) present in the msh file !" << '\n'
+            << "Only creates vega format file for the volumetric meshes.. Exiting";
+        return nullptr;
+    }
+
+    // initialize array to store the data
+    N1Delems_tag.resize(N1Delems);
+    N2Delems_tag.resize(N2Delems);
+    N3Delems_tag.resize(N3Delems);
+
+    std::string tmp_str;
+    size_t elem_tag;
+    size_t elem_type;
+    size_t elem1D_counter = 0;
+    size_t elem2D_counter = 0;
+    size_t elem3D_counter = 0;
+    std::array<size_t, 2> tmp_1arr;
+    std::array<size_t, 3> tmp_2arr;
+    std::array<size_t, 4> tmp_3arr;
+
+    while (getline(msh_stream, msh_line))
+    {
+        msh_ss.str(std::string());
+        msh_ss.clear();
+
+        if (!msh_line.empty())
+        {
+            msh_ss << msh_line;
+            msh_ss >> key_word;
+            if (!key_word.compare("$ENDELM") || !key_word.compare("$EndElements"))
+            {
+                break;
+            }
+            else
+            {
+                msh_ss.str(std::string());
+                msh_ss.clear();
+                msh_ss << msh_line;
+                msh_ss >> tmp_str;
+
+                elem_tag = stoul(tmp_str);
+                tmp_str.clear();
+                msh_ss >> tmp_str;
+                elem_type = stoul(tmp_str);
+
+                msh_ss.str(std::string());
+                msh_ss.clear();
+
+                reverse(msh_line.begin(), msh_line.end());
+                msh_ss << msh_line;
+
+                switch (elem_type)
+                {
+                    case 1: // for lines
+                        N1Delems_tag[elem1D_counter] = elem_tag;
+                        for (int jj = 1; jj >= 0; --jj)
+                        {
+                            tmp_str.clear();
+                            msh_ss >> tmp_str;
+                            reverse(tmp_str.begin(), tmp_str.end());
+                            tmp_1arr[jj] = stoul(tmp_str);
+                        }
+                        N1Delems_conn.emplace_back(tmp_1arr);
+                        ++elem1D_counter;
+                        break;
+                    case 2: // for surface elements (triangles)
+                        N2Delems_tag[elem2D_counter] = elem_tag;
+                        for (int jj = 2; jj >= 0; --jj)
+                        {
+                            tmp_str.clear();
+                            msh_ss >> tmp_str;
+                            reverse(tmp_str.begin(), tmp_str.end());
+                            tmp_2arr[jj] = stoul(tmp_str);
+                        }
+                        N2Delems_conn.emplace_back(tmp_2arr);
+                        ++elem2D_counter;
+                        break;
+                    case 4: // for volumetric elements  (tets)
+                        N3Delems_tag[elem3D_counter] = elem_tag;
+                        for (int jj = 3; jj >= 0; --jj)
+                        {
+                            tmp_str.clear();
+                            msh_ss >> tmp_str;
+                            reverse(tmp_str.begin(), tmp_str.end());
+                            tmp_3arr[jj] = stoul(tmp_str);
+                        }
+                        N3Delems_conn.emplace_back(tmp_3arr);
+                        ++elem3D_counter;
+                        break;
+                    default:
+                        break;
+                }
+            }
+        }
+    }
+    msh_stream.close();
+
+    // perform some manipulation to correct the node tags.
+    size_t tag_max = 0;
+    for (size_t ii = 0; ii < Nnodes; ++ii)
+    {
+        if (Nodes_tag[ii] > tag_max)
+        {
+            tag_max = Nodes_tag[ii];
+        }
+    }
+    tag_positions.assign(tag_max, -1);
+    for (size_t ii = 0; ii < Nnodes; ++ii)
+    {
+        tag_positions[Nodes_tag[ii] - 1] = static_cast<int>(ii); // static cast only to avoid warning
+    }
+    for (size_t ii = 0; ii < N3Delems; ++ii)
+    {
+        for (size_t jj = 0; jj < 4; ++jj)
+        {
+            N3Delems_conn[ii][jj] = tag_positions[N3Delems_conn[ii][jj] - 1];
+        }
+    }
+
+    // generate volumetric mesh
+    StdVectorOfVec3d vertices;
+    for (size_t ii = 0; ii < Nnodes; ++ii)
+    {
+        vertices.emplace_back(Nodes_xCoor[ii], Nodes_yCoor[ii], Nodes_zCoor[ii]);
+    }
+    std::vector<TetrahedralMesh::TetraArray> cells;
+    for (size_t ii = 0; ii < N3Delems; ++ii)
+    {
+        cells.emplace_back(N3Delems_conn[ii]);
+    }
+    auto tetMesh = std::make_shared<TetrahedralMesh>();
+    tetMesh->initialize(vertices, cells, false);
+    return tetMesh;
+};
+
+}
\ No newline at end of file
diff --git a/Base/Geometry/Reader/imstkMSHMeshIO.h b/Base/Geometry/Reader/imstkMSHMeshIO.h
new file mode 100644
index 0000000000000000000000000000000000000000..715b43b2ddfdafbd1931e159c4c1e189b2b78f15
--- /dev/null
+++ b/Base/Geometry/Reader/imstkMSHMeshIO.h
@@ -0,0 +1,63 @@
+/*=========================================================================
+
+    Library: iMSTK
+
+    Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+    & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0.txt
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+
+    =========================================================================*/
+
+#ifndef imstkMSHMeshIO_h
+#define imstkMSHMeshIO_h
+
+#include <memory>
+#include <vector>
+#include <array>
+
+#include "imstkMeshIO.h"
+#include "imstkVolumetricMesh.h"
+
+namespace imstk
+{
+
+///
+/// \class MSHMeshIO
+///
+/// \brief Contains utility to generate imstk::VolumetricMesh from mesh given in msh
+/// file format. Only works for tet meshes.
+///
+class MSHMeshIO
+{
+public:
+    ///
+    /// \brief Defualt Constructor
+    ///
+    MSHMeshIO() = default;
+
+    ///
+    /// \brief Default Destructor
+    ///
+    ~MSHMeshIO() = default;
+
+    ///
+    /// \brief Read and generate a volumetric mesh given a external msh file
+    ///
+    static std::shared_ptr<imstk::VolumetricMesh> read(const std::string &filePath,const MeshFileType meshType);
+
+};
+
+}
+
+#endif// imstkMSHMeshIO_h
\ No newline at end of file
diff --git a/Base/Geometry/Reader/imstkMeshReader.cpp b/Base/Geometry/Reader/imstkMeshIO.cpp
similarity index 63%
rename from Base/Geometry/Reader/imstkMeshReader.cpp
rename to Base/Geometry/Reader/imstkMeshIO.cpp
index 6035660ee7e00b7d01a86786bb3e4dd2da575b81..8a1ae90d1fee27cacaf31c8d833bf56aa86a09ae 100644
--- a/Base/Geometry/Reader/imstkMeshReader.cpp
+++ b/Base/Geometry/Reader/imstkMeshIO.cpp
@@ -19,12 +19,12 @@
 
    =========================================================================*/
 
-#include "imstkMeshReader.h"
-
 #include <sys/stat.h>
 
-#include "imstkVTKMeshReader.h"
-#include "imstkVegaMeshReader.h"
+#include "imstkMeshIO.h"
+#include "imstkVTKMeshIO.h"
+#include "imstkVegaMeshIO.h"
+#include "imstkMSHMeshIO.h"
 
 #include "g3log/g3log.hpp"
 
@@ -32,15 +32,15 @@ namespace imstk
 {
 
 std::shared_ptr<Mesh>
-MeshReader::read(const std::string& filePath)
+MeshIO::read(const std::string& filePath)
 {
-    if (!MeshReader::fileExists(filePath))
+    if (!MeshIO::fileExists(filePath))
     {
-        LOG(WARNING) << "MeshReader::read error: file not found: " << filePath;
+        LOG(WARNING) << "MeshIO::read error: file not found: " << filePath;
         return nullptr;
     }
 
-    MeshFileType meshType = MeshReader::getFileType(filePath);
+    MeshFileType meshType = MeshIO::getFileType(filePath);
     switch (meshType)
     {
     case MeshFileType::VTK :
@@ -49,33 +49,36 @@ MeshReader::read(const std::string& filePath)
     case MeshFileType::STL :
     case MeshFileType::PLY :
     case MeshFileType::OBJ :
-        return VTKMeshReader::read(filePath, meshType);
+        return VTKMeshIO::read(filePath, meshType);
         break;
     case MeshFileType::VEG :
-        return VegaMeshReader::read(filePath, meshType);
+        return VegaMeshIO::read(filePath, meshType);
+        break;
+    case MeshFileType::MSH:
+        return MSHMeshIO::read(filePath, meshType);
         break;
     }
 
-    LOG(WARNING) << "MeshReader::read error: file type not supported";
+    LOG(WARNING) << "MeshIO::read error: file type not supported";
     return nullptr;
 }
 
 bool
-MeshReader::fileExists(const std::string& file)
+MeshIO::fileExists(const std::string& file)
 {
     struct stat buf;
     return (stat(file.c_str(), &buf) == 0);
 }
 
 const MeshFileType
-MeshReader::getFileType(const std::string& filePath)
+MeshIO::getFileType(const std::string& filePath)
 {
     MeshFileType meshType = MeshFileType::UNKNOWN;
 
     std::string extString = filePath.substr(filePath.find_last_of(".") + 1);
     if (extString.empty())
     {
-        LOG(WARNING) << "MeshReader::getFileType error: invalid file name";
+        LOG(WARNING) << "MeshIO::getFileType error: invalid file name";
         return meshType;
     }
 
@@ -107,12 +110,31 @@ MeshReader::getFileType(const std::string& filePath)
     {
         meshType = MeshFileType::VEG;
     }
+    else if (extString == "msh" || extString == "MSH")
+    {
+        meshType = MeshFileType::MSH;
+    }
     else
     {
-        LOG(WARNING) << "MeshReader::getFileType error: unknown file extension";
+        LOG(WARNING) << "MeshIO::getFileType error: unknown file extension";
     }
 
     return meshType;
 }
 
+bool
+MeshIO::write(const std::shared_ptr<imstk::Mesh> imstkMesh, const std::string& filePath)
+{
+    MeshFileType meshType = MeshIO::getFileType(filePath);
+    switch (meshType)
+    {
+        case MeshFileType::VEG:
+            return VegaMeshIO::write(imstkMesh, filePath, meshType);
+            break;
+    }
+
+    LOG(WARNING) << "MeshIO::write error: file type not supported";
+    return nullptr;
+}
+
 } // imstk
diff --git a/Base/Geometry/Reader/imstkMeshReader.h b/Base/Geometry/Reader/imstkMeshIO.h
similarity index 82%
rename from Base/Geometry/Reader/imstkMeshReader.h
rename to Base/Geometry/Reader/imstkMeshIO.h
index dc65839cd35d3e3ac5dc27c66a7a68d9f61a69bf..d6aafbf3625e939f4a20e98423129eb71f1fe8d9 100644
--- a/Base/Geometry/Reader/imstkMeshReader.h
+++ b/Base/Geometry/Reader/imstkMeshIO.h
@@ -19,8 +19,8 @@
 
    =========================================================================*/
 
-#ifndef imstkMeshReader_h
-#define imstkMeshReader_h
+#ifndef imstkMeshIO_h
+#define imstkMeshIO_h
 
 // std library
 #include <memory>
@@ -45,39 +45,46 @@ enum MeshFileType
     STL,
     PLY,
     OBJ,
-    VEG
+    VEG,
+    MSH
 };
 
 ///
-/// \class MeshReader
+/// \class MeshIO
 ///
-/// \brief Mesh data reader
+/// \brief Mesh data IO
 ///
-class MeshReader
+class MeshIO
 {
 public:
 
     ///
     /// \brief Constructor
     ///
-    MeshReader() = default;
+    MeshIO() = default;
 
     ///
     /// \brief Destructor
     ///
-    ~MeshReader() = default;
+    ~MeshIO() = default;
 
     ///
     /// \brief Read external file
     ///
     static std::shared_ptr<Mesh> read(const std::string& filePath);
 
+    ///
+    /// \brief Write external file
+    ///
+    static bool write(const std::shared_ptr<imstk::Mesh> imstkMesh, const std::string& filePath);
+
     ///
     /// \brief Returns true if the file exists, else false
     ///
     static bool fileExists(const std::string& file);
 
 protected:
+
     ///
     /// \brief Returns the type of the file
     ///
@@ -87,4 +94,4 @@ protected:
 
 } // imstk
 
-#endif // ifndef imstkMeshReader_h
+#endif // ifndef imstkMeshIO_h
diff --git a/Base/Geometry/Reader/imstkVTKMeshReader.cpp b/Base/Geometry/Reader/imstkVTKMeshIO.cpp
similarity index 67%
rename from Base/Geometry/Reader/imstkVTKMeshReader.cpp
rename to Base/Geometry/Reader/imstkVTKMeshIO.cpp
index c543664b66a0a86c5b85b2b0d29e39e9260803f1..0b99f68886476d5db7bdd17e7f887378e57a2536 100644
--- a/Base/Geometry/Reader/imstkVTKMeshReader.cpp
+++ b/Base/Geometry/Reader/imstkVTKMeshIO.cpp
@@ -19,7 +19,7 @@
 
    =========================================================================*/
 
-#include "imstkVTKMeshReader.h"
+#include "imstkVTKMeshIO.h"
 
 #include "vtkSmartPointer.h"
 #include "vtkGenericDataObjectReader.h"
@@ -36,43 +36,43 @@ namespace imstk
 {
 
 std::shared_ptr<Mesh>
-VTKMeshReader::read(const std::string& filePath, MeshFileType meshType)
+VTKMeshIO::read(const std::string& filePath, MeshFileType meshType)
 {
     switch (meshType)
     {
     case MeshFileType::VTK :
     {
-        return VTKMeshReader::readVtkGenericFormatData<vtkGenericDataObjectReader>(filePath);
+        return VTKMeshIO::readVtkGenericFormatData<vtkGenericDataObjectReader>(filePath);
         break;
     }
     case MeshFileType::VTU :
     {
-        return VTKMeshReader::readVtkUnstructuredGrid<vtkXMLUnstructuredGridReader>(filePath);
+        return VTKMeshIO::readVtkUnstructuredGrid<vtkXMLUnstructuredGridReader>(filePath);
         break;
     }
     case MeshFileType::VTP :
     {
-        return VTKMeshReader::readVtkPolyData<vtkXMLPolyDataReader>(filePath);
+        return VTKMeshIO::readVtkPolyData<vtkXMLPolyDataReader>(filePath);
         break;
     }
     case MeshFileType::STL :
     {
-        return VTKMeshReader::readVtkPolyData<vtkSTLReader>(filePath);
+        return VTKMeshIO::readVtkPolyData<vtkSTLReader>(filePath);
         break;
     }
     case MeshFileType::PLY :
     {
-        return VTKMeshReader::readVtkPolyData<vtkPLYReader>(filePath);
+        return VTKMeshIO::readVtkPolyData<vtkPLYReader>(filePath);
         break;
     }
     case MeshFileType::OBJ :
     {
-        return VTKMeshReader::readVtkPolyData<vtkOBJReader>(filePath);
+        return VTKMeshIO::readVtkPolyData<vtkOBJReader>(filePath);
         break;
     }
     default :
     {
-        LOG(WARNING) << "VTKMeshReader::read error: file type not supported";
+        LOG(WARNING) << "VTKMeshIO::read error: file type not supported";
         break;
     }
     }
@@ -80,7 +80,7 @@ VTKMeshReader::read(const std::string& filePath, MeshFileType meshType)
 
 template<typename ReaderType>
 std::shared_ptr<Mesh>
-VTKMeshReader::readVtkGenericFormatData(const std::string& filePath)
+VTKMeshIO::readVtkGenericFormatData(const std::string& filePath)
 {
     auto reader = vtkSmartPointer<ReaderType>::New();
     reader->SetFileName(filePath.c_str());
@@ -88,64 +88,64 @@ VTKMeshReader::readVtkGenericFormatData(const std::string& filePath)
 
     if (vtkPolyData* vtkMesh = reader->GetPolyDataOutput())
     {
-        return VTKMeshReader::convertVtkPolyDataToSurfaceMesh(vtkMesh);
+        return VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkMesh);
     }
     else if (vtkUnstructuredGrid* vtkMesh = reader->GetUnstructuredGridOutput())
     {
-        return VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
+        return VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
     }
     else
     {
-        LOG(WARNING) << "VTKMeshReader::readVtkGenericFormatData error: could not read with VTK reader.";
+        LOG(WARNING) << "VTKMeshIO::readVtkGenericFormatData error: could not read with VTK reader.";
         return nullptr;
     }
 }
 
 template<typename ReaderType>
 std::shared_ptr<SurfaceMesh>
-VTKMeshReader::readVtkPolyData(const std::string& filePath)
+VTKMeshIO::readVtkPolyData(const std::string& filePath)
 {
     auto reader = vtkSmartPointer<ReaderType>::New();
     reader->SetFileName(filePath.c_str());
     reader->Update();
 
     vtkPolyData* vtkMesh = reader->GetOutput();
-    return VTKMeshReader::convertVtkPolyDataToSurfaceMesh(vtkMesh);
+    return VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkMesh);
 }
 
 template<typename ReaderType>
 std::shared_ptr<VolumetricMesh>
-VTKMeshReader::readVtkUnstructuredGrid(const std::string& filePath)
+VTKMeshIO::readVtkUnstructuredGrid(const std::string& filePath)
 {
     auto reader = vtkSmartPointer<ReaderType>::New();
     reader->SetFileName(filePath.c_str());
     reader->Update();
 
     vtkUnstructuredGrid* vtkMesh = reader->GetOutput();
-    return VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
+    return VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkMesh);
 }
 
 std::shared_ptr<SurfaceMesh>
-VTKMeshReader::convertVtkPolyDataToSurfaceMesh(vtkPolyData* vtkMesh)
+VTKMeshIO::convertVtkPolyDataToSurfaceMesh(vtkPolyData* vtkMesh)
 {
     if(!vtkMesh)
     {
-        LOG(WARNING) << "VTKMeshReader::convertVtkPolyDataToSurfaceMesh error: could not read with VTK reader.";
+        LOG(WARNING) << "VTKMeshIO::convertVtkPolyDataToSurfaceMesh error: could not read with VTK reader.";
         return nullptr;
     }
 
     StdVectorOfVec3d vertices;
-    VTKMeshReader::copyVertices(vtkMesh->GetPoints(), vertices);
+    VTKMeshIO::copyVertices(vtkMesh->GetPoints(), vertices);
 
     std::vector<SurfaceMesh::TriangleArray> triangles;
-    VTKMeshReader::copyCells<3>(vtkMesh->GetPolys(), triangles);
+    VTKMeshIO::copyCells<3>(vtkMesh->GetPolys(), triangles);
 
     auto mesh = std::make_shared<SurfaceMesh>();
     mesh->initialize(vertices, triangles, true);
 
     // Point Data
     std::map<std::string, StdVectorOfVectorf> dataMap;
-    VTKMeshReader::copyPointData(vtkMesh->GetPointData(), dataMap);
+    VTKMeshIO::copyPointData(vtkMesh->GetPointData(), dataMap);
     if (!dataMap.empty())
     {
       mesh->setPointDataMap(dataMap);
@@ -162,22 +162,22 @@ VTKMeshReader::convertVtkPolyDataToSurfaceMesh(vtkPolyData* vtkMesh)
 }
 
 std::shared_ptr<VolumetricMesh>
-VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* vtkMesh)
+VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* vtkMesh)
 {
     if(!vtkMesh)
     {
-        LOG(WARNING) << "VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh error: could not read with VTK reader.";
+        LOG(WARNING) << "VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh error: could not read with VTK reader.";
         return nullptr;
     }
 
     StdVectorOfVec3d vertices;
-    VTKMeshReader::copyVertices(vtkMesh->GetPoints(), vertices);
+    VTKMeshIO::copyVertices(vtkMesh->GetPoints(), vertices);
 
     int cellType = vtkMesh->GetCellType(vtkMesh->GetNumberOfCells()-1);
     if( cellType == VTK_TETRA )
     {
         std::vector<TetrahedralMesh::TetraArray> cells;
-        VTKMeshReader::copyCells<4>(vtkMesh->GetCells(), cells);
+        VTKMeshIO::copyCells<4>(vtkMesh->GetCells(), cells);
 
         auto mesh = std::make_shared<TetrahedralMesh>();
         mesh->initialize(vertices, cells, false);
@@ -186,7 +186,7 @@ VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* v
     else if( cellType == VTK_HEXAHEDRON )
     {
         std::vector<HexahedralMesh::HexaArray> cells;
-        VTKMeshReader::copyCells<8>(vtkMesh->GetCells(), cells);
+        VTKMeshIO::copyCells<8>(vtkMesh->GetCells(), cells);
 
         auto mesh = std::make_shared<HexahedralMesh>();
         mesh->initialize(vertices, cells, false);
@@ -194,18 +194,18 @@ VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* v
     }
     else
     {
-        LOG(WARNING) << "VTKMeshReader::convertVtkUnstructuredGridToVolumetricMesh error: No support for vtkCellType="
+        LOG(WARNING) << "VTKMeshIO::convertVtkUnstructuredGridToVolumetricMesh error: No support for vtkCellType="
                      << cellType << ".";
         return nullptr;
     }
 }
 
 void
-VTKMeshReader::copyVertices(vtkPoints* points, StdVectorOfVec3d& vertices)
+VTKMeshIO::copyVertices(vtkPoints* points, StdVectorOfVec3d& vertices)
 {
     if(!points)
     {
-        LOG(WARNING) << "VTKMeshReader::copyVertices error: No points found.";
+        LOG(WARNING) << "VTKMeshIO::copyVertices error: No points found.";
         return;
     }
 
@@ -219,11 +219,11 @@ VTKMeshReader::copyVertices(vtkPoints* points, StdVectorOfVec3d& vertices)
 
 template<size_t dim>
 void
-VTKMeshReader::copyCells(vtkCellArray* vtkCells, std::vector<std::array<size_t,dim>>& cells)
+VTKMeshIO::copyCells(vtkCellArray* vtkCells, std::vector<std::array<size_t,dim>>& cells)
 {
     if(!vtkCells)
     {
-        LOG(WARNING) << "VTKMeshReader::copyCells error: No cells found.";
+        LOG(WARNING) << "VTKMeshIO::copyCells error: No cells found.";
         return;
     }
 
@@ -245,7 +245,7 @@ VTKMeshReader::copyCells(vtkCellArray* vtkCells, std::vector<std::array<size_t,d
 }
 
 void
-VTKMeshReader::copyPointData(vtkPointData* pointData, std::map<std::string, StdVectorOfVectorf>& dataMap)
+VTKMeshIO::copyPointData(vtkPointData* pointData, std::map<std::string, StdVectorOfVectorf>& dataMap)
 {
     if(!pointData)
     {
diff --git a/Base/Geometry/Reader/imstkVTKMeshReader.h b/Base/Geometry/Reader/imstkVTKMeshIO.h
similarity index 91%
rename from Base/Geometry/Reader/imstkVTKMeshReader.h
rename to Base/Geometry/Reader/imstkVTKMeshIO.h
index 24e0b45b3246c3aeb0f2c0c871c3555c9a199beb..0ca2f6d7d4a67425def246417c534a49d9fb7a70 100644
--- a/Base/Geometry/Reader/imstkVTKMeshReader.h
+++ b/Base/Geometry/Reader/imstkVTKMeshIO.h
@@ -19,8 +19,8 @@
 
    =========================================================================*/
 
-#ifndef imstkVTKMeshReader_h
-#define imstkVTKMeshReader_h
+#ifndef imstkVTKMeshIO_h
+#define imstkVTKMeshIO_h
 
 #include <memory>
 
@@ -30,7 +30,7 @@
 #include "vtkCellArray.h"
 #include "vtkPointData.h"
 
-#include "imstkMeshReader.h"
+#include "imstkMeshIO.h"
 #include "imstkSurfaceMesh.h"
 #include "imstkTetrahedralMesh.h"
 #include "imstkHexahedralMesh.h"
@@ -39,23 +39,23 @@ namespace imstk
 {
 
 ///
-/// \class VTKMeshReader
+/// \class VTKMeshIO
 ///
 /// \brief
 ///
-class VTKMeshReader
+class VTKMeshIO
 {
 public:
 
     ///
     /// \brief Default constructor
     ///
-    VTKMeshReader() = default;
+    VTKMeshIO() = default;
 
     ///
     /// \brief Default destructor
     ///
-    ~VTKMeshReader() = default;
+    ~VTKMeshIO() = default;
 
     ///
     /// \brief
@@ -112,4 +112,4 @@ protected:
 
 } // imstk
 
-#endif // ifndef imstkVTKMeshReader_h
+#endif // ifndef imstkVTKMeshIO_h
diff --git a/Base/Geometry/Reader/imstkVegaMeshIO.cpp b/Base/Geometry/Reader/imstkVegaMeshIO.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..753a53154002367967f01f31d0ff826d02032482
--- /dev/null
+++ b/Base/Geometry/Reader/imstkVegaMeshIO.cpp
@@ -0,0 +1,216 @@
+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+   =========================================================================*/
+
+#include "imstkVegaMeshIO.h"
+#include "imstkMeshIO.h"
+#include "imstkTetrahedralMesh.h"
+#include "imstkHexahedralMesh.h"
+
+#include "g3log/g3log.hpp"
+
+#include "tetMesh.h"
+
+namespace imstk
+{
+
+std::shared_ptr<VolumetricMesh>
+VegaMeshIO::read(const std::string& filePath, MeshFileType meshType)
+{
+    if (meshType != MeshFileType::VEG)
+    {
+        LOG(WARNING) << "VegaMeshIO::read error: file type not supported";
+        return nullptr;
+    }
+
+    // Read Vega Mesh
+    std::shared_ptr<vega::VolumetricMesh> vegaMesh = VegaMeshIO::readVegaMesh(filePath);
+
+    // Convert to Volumetric Mesh
+    return VegaMeshIO::convertVegaMeshToVolumetricMesh(vegaMesh);
+}
+
+bool
+VegaMeshIO::write(const std::shared_ptr<imstk::Mesh> imstkMesh, const std::string& filePath, const MeshFileType meshType)
+{
+    if (meshType != MeshFileType::VEG)
+    {
+        LOG(WARNING) << "VegaMeshIO::write error: file type is not veg";
+        return false;
+    }
+    // extract volumetric mesh
+    const auto imstkVolMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(imstkMesh);
+    if (!imstkVolMesh)
+    {
+        LOG(WARNING) << "VegaMeshIO::write error: imstk::Mesh is not a volumetric mesh";
+        return false;
+    }
+
+    switch (imstkVolMesh->getType())
+    {
+        case Geometry::Type::TetrahedralMesh:
+        case Geometry::Type::HexahedralMesh:
+            auto vegaMesh = convertVolumetricMeshToVegaMesh(imstkVolMesh);
+            if (vegaMesh==nullptr)
+            {
+                LOG(WARNING) << "VegaMeshIO::write error: failed to convert volumetric mesh to vega mesh";
+                return false;
+            }
+            const auto fileName = const_cast<char *>(filePath.c_str());
+            const int write_status = vegaMesh->save(fileName);
+            if (write_status != 0)
+            {
+                LOG(WARNING) << "VegaMeshIO::write error: failed to write .veg file";
+                return false;
+            }
+            return true;
+            break;
+    }
+
+    LOG(WARNING) << "VegaMeshIO::write error: this element type not supported in vega";
+    return false;
+}
+
+std::shared_ptr<vega::VolumetricMesh>
+VegaMeshIO::readVegaMesh(const std::string& filePath)
+{
+    auto fileName = const_cast<char*>(filePath.c_str());
+    std::shared_ptr<vega::VolumetricMesh> vegaMesh(vega::VolumetricMeshLoader::load(fileName));
+    return vegaMesh;
+}
+
+std::shared_ptr<imstk::VolumetricMesh>
+VegaMeshIO::convertVegaMeshToVolumetricMesh(std::shared_ptr<vega::VolumetricMesh> vegaMesh)
+{
+    // Copy vertices
+    StdVectorOfVec3d vertices;
+    VegaMeshIO::copyVertices(vegaMesh, vertices);
+
+    // Copy cells
+    auto cellType = vegaMesh->getElementType();
+    std::shared_ptr<imstk::VolumetricMesh> mesh;
+    if(cellType == vega::VolumetricMesh::TET)
+    {
+        std::vector<TetrahedralMesh::TetraArray> cells;
+        VegaMeshIO::copyCells<4>(vegaMesh, cells);
+
+        auto tetMesh = std::make_shared<TetrahedralMesh>();
+        tetMesh->initialize(vertices, cells, false);
+        mesh = tetMesh;
+    }
+    else if(cellType == vega::VolumetricMesh::CUBIC)
+    {
+        std::vector<HexahedralMesh::HexaArray> cells;
+        VegaMeshIO::copyCells<8>(vegaMesh, cells);
+
+        auto hexMesh = std::make_shared<HexahedralMesh>();
+        hexMesh->initialize(vertices, cells, false);
+        mesh = hexMesh;
+    }
+    else
+    {
+        vegaMesh.reset();
+        LOG(WARNING) << "VegaMeshIO::read error: invalid cell type.";
+        return nullptr;
+    }
+
+    // Keep track of the vega mesh to initialize dynamical model
+    mesh->setAttachedVegaMesh(vegaMesh);
+    return mesh;
+}
+
+void
+VegaMeshIO::copyVertices(std::shared_ptr<vega::VolumetricMesh> vegaMesh,
+                             StdVectorOfVec3d& vertices)
+{
+    for(size_t i = 0; i < vegaMesh->getNumVertices(); ++i)
+    {
+        auto pos = *vegaMesh->getVertex(i);
+        vertices.emplace_back(pos[0], pos[1], pos[2]);
+    }
+}
+
+template<size_t dim>
+void
+VegaMeshIO::copyCells(std::shared_ptr<vega::VolumetricMesh> vegaMesh,
+                          std::vector<std::array<size_t,dim>>& cells)
+{
+    std::array<size_t,dim> cell;
+    for(size_t cellId = 0; cellId < vegaMesh->getNumElements(); ++cellId)
+    {
+        for (size_t i = 0; i < vegaMesh->getNumElementVertices(); ++i)
+        {
+            cell[i] = vegaMesh->getVertexIndex(cellId,i);
+        }
+        cells.emplace_back(cell);
+    }
+}
+
+std::shared_ptr<vega::VolumetricMesh>
+VegaMeshIO::convertVolumetricMeshToVegaMesh(const std::shared_ptr<imstk::VolumetricMesh> imstkVolMesh)
+{
+     // as of now, only works for TET elements
+	if (imstkVolMesh->getType() == Geometry::Type::TetrahedralMesh)
+    {
+        // Using default material properties to append to the .veg file
+        const double E=1E6;
+        const double nu=0.45;
+        const double density=1000.0;
+
+        auto imstkVolTetMesh = std::dynamic_pointer_cast<imstk::TetrahedralMesh>(imstkVolMesh);
+
+        auto vertexArray = imstkVolMesh->getVertexPositions();
+        std::vector<double> vertices;
+        for (const auto & node : vertexArray)
+        {
+            vertices.emplace_back(node(0));
+            vertices.emplace_back(node(1));
+            vertices.emplace_back(node(2));
+        }
+
+        auto tetArray = imstkVolTetMesh->getTetrahedraVertices();
+        std::vector<int> elements;
+        for (const auto & tet : tetArray)
+        {
+            elements.emplace_back(tet[0]);
+            elements.emplace_back(tet[1]);
+            elements.emplace_back(tet[2]);
+            elements.emplace_back(tet[3]);
+        }
+
+        std::shared_ptr<vega::VolumetricMesh> vegaMesh(new vega::TetMesh(imstkVolTetMesh->getNumVertices(), &vertices[0], imstkVolTetMesh->getNumTetrahedra(), &elements[0], E, nu, density));
+        if (!vegaMesh)
+        {
+            LOG(WARNING) << "VegaMeshIO::convertVolumetricMeshToVegaMesh error: Failed to create vega mesh";
+            return nullptr;
+        }
+        else
+        {
+            return vegaMesh;
+        }
+    }
+    else
+    {
+        LOG(WARNING) << "VegaMeshIO::convertVolumetricMeshToVegaMesh error: Geometry type not supported";
+        return nullptr;
+    }
+}
+
+} // imstk
diff --git a/Base/Geometry/Reader/imstkVegaMeshReader.h b/Base/Geometry/Reader/imstkVegaMeshIO.h
similarity index 82%
rename from Base/Geometry/Reader/imstkVegaMeshReader.h
rename to Base/Geometry/Reader/imstkVegaMeshIO.h
index 25f950c44a08a3210d21c4773232beebba5539dc..2cac5a1e9ffe0a3096221e360dec29d242f15266 100644
--- a/Base/Geometry/Reader/imstkVegaMeshReader.h
+++ b/Base/Geometry/Reader/imstkVegaMeshIO.h
@@ -19,12 +19,12 @@
 
    =========================================================================*/
 
-#ifndef imstkVegaMeshReader_h
-#define imstkVegaMeshReader_h
+#ifndef imstkVegaMeshIO_h
+#define imstkVegaMeshIO_h
 
 #include <memory>
 
-#include "imstkMeshReader.h"
+#include "imstkMeshIO.h"
 #include "imstkVolumetricMesh.h"
 
 // Vega
@@ -35,23 +35,24 @@ namespace imstk
 {
 
 ///
-/// \class VegaMeshReader
+/// \class VegaMeshIO
 ///
 /// \brief Contains utility classes that convert vega volume mesh to volume mesh and
 /// vice-versa
 ///
-class VegaMeshReader
+class VegaMeshIO
 {
+
 public:
     ///
     /// \brief Default constructor
     ///
-    VegaMeshReader() = default;
+    VegaMeshIO() = default;
 
     ///
     /// \brief Default destructor
     ///
-    ~VegaMeshReader() = default;
+    ~VegaMeshIO() = default;
 
     ///
     /// \brief Read and generate volumetric mesh given a external vega mesh file
@@ -63,6 +64,11 @@ public:
     ///
     static std::shared_ptr<vega::VolumetricMesh> readVegaMesh(const std::string& filePath);
 
+    ///
+    /// \brief Write a volumetric mesh in vega file format
+    ///
+    static bool write(const std::shared_ptr<imstk::Mesh> imstkMesh, const std::string& filePath, const MeshFileType meshType);
+
 protected:
     ///
     /// \brief Generate volumetric mesh given a vega volume mesh
@@ -72,7 +78,7 @@ protected:
     ///
     /// \brief Generate a vega volume mesh given volumetric mesh
     ///
-    static std::shared_ptr<vega::VolumetricMesh> convertVolumetricMeshToVegaMesh(std::shared_ptr<VolumetricMesh> volumeMesh);
+    static std::shared_ptr<vega::VolumetricMesh> convertVolumetricMeshToVegaMesh(const std::shared_ptr<imstk::VolumetricMesh> volumeMesh);
 
     ///
     /// \brief
@@ -88,4 +94,4 @@ protected:
 
 } // imstk
 
-#endif // ifndef imstkVegaMeshReader_h
+#endif // ifndef imstkVegaMeshIO_h
diff --git a/Base/Geometry/Reader/imstkVegaMeshReader.cpp b/Base/Geometry/Reader/imstkVegaMeshReader.cpp
deleted file mode 100644
index bd621fcf2aba3c64e09bd539cea538e9eb6ec722..0000000000000000000000000000000000000000
--- a/Base/Geometry/Reader/imstkVegaMeshReader.cpp
+++ /dev/null
@@ -1,124 +0,0 @@
-/*=========================================================================
-
-   Library: iMSTK
-
-   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
-   & Imaging in Medicine, Rensselaer Polytechnic Institute.
-
-   Licensed under the Apache License, Version 2.0 (the "License");
-   you may not use this file except in compliance with the License.
-   You may obtain a copy of the License at
-
-      http://www.apache.org/licenses/LICENSE-2.0.txt
-
-   Unless required by applicable law or agreed to in writing, software
-   distributed under the License is distributed on an "AS IS" BASIS,
-   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-   See the License for the specific language governing permissions and
-   limitations under the License.
-
-   =========================================================================*/
-
-#include "imstkVegaMeshReader.h"
-#include "imstkMeshReader.h"
-#include "imstkTetrahedralMesh.h"
-#include "imstkHexahedralMesh.h"
-
-#include "g3log/g3log.hpp"
-
-namespace imstk
-{
-
-std::shared_ptr<VolumetricMesh>
-VegaMeshReader::read(const std::string& filePath, MeshFileType meshType)
-{
-    if (meshType != MeshFileType::VEG)
-    {
-        LOG(WARNING) << "VegaMeshReader::read error: file type not supported";
-        return nullptr;
-    }
-
-    // Read Vega Mesh
-    std::shared_ptr<vega::VolumetricMesh> vegaMesh = VegaMeshReader::readVegaMesh(filePath);
-
-    // Convert to Volumetric Mesh
-    return VegaMeshReader::convertVegaMeshToVolumetricMesh(vegaMesh);
-}
-
-
-std::shared_ptr<vega::VolumetricMesh>
-VegaMeshReader::readVegaMesh(const std::string& filePath)
-{
-    auto fileName = const_cast<char*>(filePath.c_str());
-    std::shared_ptr<vega::VolumetricMesh> vegaMesh(vega::VolumetricMeshLoader::load(fileName));
-    return vegaMesh;
-}
-
-std::shared_ptr<imstk::VolumetricMesh>
-VegaMeshReader::convertVegaMeshToVolumetricMesh(std::shared_ptr<vega::VolumetricMesh> vegaMesh)
-{
-    // Copy vertices
-    StdVectorOfVec3d vertices;
-    VegaMeshReader::copyVertices(vegaMesh, vertices);
-
-    // Copy cells
-    auto cellType = vegaMesh->getElementType();
-    std::shared_ptr<imstk::VolumetricMesh> mesh;
-    if(cellType == vega::VolumetricMesh::TET)
-    {
-        std::vector<TetrahedralMesh::TetraArray> cells;
-        VegaMeshReader::copyCells<4>(vegaMesh, cells);
-
-        auto tetMesh = std::make_shared<TetrahedralMesh>();
-        tetMesh->initialize(vertices, cells, false);
-        mesh = tetMesh;
-    }
-    else if(cellType == vega::VolumetricMesh::CUBIC)
-    {
-        std::vector<HexahedralMesh::HexaArray> cells;
-        VegaMeshReader::copyCells<8>(vegaMesh, cells);
-
-        auto hexMesh = std::make_shared<HexahedralMesh>();
-        hexMesh->initialize(vertices, cells, false);
-        mesh = hexMesh;
-    }
-    else
-    {
-        vegaMesh.reset();
-        LOG(WARNING) << "VegaMeshReader::read error: invalid cell type.";
-        return nullptr;
-    }
-
-    // Keep track of the vega mesh to initialize dynamical model
-    mesh->setAttachedVegaMesh(vegaMesh);
-    return mesh;
-}
-
-void
-VegaMeshReader::copyVertices(std::shared_ptr<vega::VolumetricMesh> vegaMesh,
-                             StdVectorOfVec3d& vertices)
-{
-    for(size_t i = 0; i < vegaMesh->getNumVertices(); ++i)
-    {
-        auto pos = *vegaMesh->getVertex(i);
-        vertices.emplace_back(pos[0], pos[1], pos[2]);
-    }
-}
-
-template<size_t dim>
-void
-VegaMeshReader::copyCells(std::shared_ptr<vega::VolumetricMesh> vegaMesh,
-                          std::vector<std::array<size_t,dim>>& cells)
-{
-    std::array<size_t,dim> cell;
-    for(size_t cellId = 0; cellId < vegaMesh->getNumElements(); ++cellId)
-    {
-        for (size_t i = 0; i < vegaMesh->getNumElementVertices(); ++i)
-        {
-            cell[i] = vegaMesh->getVertexIndex(cellId,i);
-        }
-        cells.emplace_back(cell);
-    }
-}
-
-} // imstk
diff --git a/Examples/Sandbox/main.cpp b/Examples/Sandbox/main.cpp
index f0a5ea97311cf15f8adf06cdac91a7cb9226aec8..cb5ea6dc61562da44e6c0f12b3d1b605bb7f3ea0 100644
--- a/Examples/Sandbox/main.cpp
+++ b/Examples/Sandbox/main.cpp
@@ -1,4 +1,4 @@
-#define DATA_ROOT_PATH "." // Change to your data ressource directory
+#define DATA_ROOT_PATH "M:/Mohit Tyagi/iMSTK/build/Innerbuild/Examples/Sandbox/Resources" // Change to your data ressource directory
 
 #include <cstring>
 #include <iostream>
@@ -33,7 +33,7 @@
 #include "imstkCube.h"
 #include "imstkTetrahedralMesh.h"
 #include "imstkSurfaceMesh.h"
-#include "imstkMeshReader.h"
+#include "imstkMeshIO.h"
 #include "imstkLineMesh.h"
 
 // Maps
@@ -96,6 +96,7 @@ void testPbdVolume();
 void testPbdCloth();
 void testPbdCollision();
 void testLineMesh();
+void testMshAndVegaIO();
 
 int main()
 {
@@ -124,13 +125,74 @@ int main()
     //testTwoOmnis();
     //testVectorPlotters();
     //testPbdVolume();
-    testPbdCloth();
+    //testPbdCloth();
     //testPbdCollision();
     //testLineMesh();
+    testMshAndVegaIO();
 
     return 0;
 }
 
+void testMshAndVegaIO()
+{
+    // SDK and Scene
+    auto sdk = std::make_shared<imstk::SimulationManager>();
+    auto scene = sdk->createNewScene("SceneTestMesh");
+
+    // Load a volumetric mesh (from .msh file)
+    std::string ifile = DATA_ROOT_PATH"/sofa_msh_files/liver2.msh";
+    auto volMeshA = imstk::MeshIO::read(ifile);
+    if (!volMeshA)
+    {
+        LOG(WARNING) << "Failed to read msh file : " << ifile;
+        return;
+    }
+
+    // Extract surface mesh
+    auto volumeMeshA = std::dynamic_pointer_cast<imstk::VolumetricMesh>(volMeshA); // change to any volumetric mesh above
+    volumeMeshA->computeAttachedSurfaceMesh();
+    auto surfaceMeshA = volumeMeshA->getAttachedSurfaceMesh();
+
+    // Create object A
+    auto objectA = std::make_shared<imstk::VisualObject>("meshObjectMSH");
+    objectA->setVisualGeometry(surfaceMeshA);
+
+    // Write a .veg file
+    std::string ofile = DATA_ROOT_PATH"/VEGA_veg_files/liver2.veg";
+    auto writeStatus = imstk::MeshIO::write(volMeshA, ofile);
+    std::cout << "------------------------------Summary----------------------------------------------------\n";
+    std::cout << "Following file converion: " << ((writeStatus) ? "Success \n" : "Failure \n");
+    std::cout << "\n Input mesh file : \n" << ifile << std::endl;
+    std::cout << "\n Output mesh file: \n" << ofile << std::endl;
+
+    // Read the above written veg file
+    auto volMeshB = imstk::MeshIO::read(ofile);
+    if (!volMeshB)
+    {
+        LOG(WARNING) << "Failed to extract topology/geometry from the veg file : " << ofile;
+        return;
+    }
+
+    // Extract surface mesh
+    auto volumeMeshB = std::dynamic_pointer_cast<imstk::VolumetricMesh>(volMeshB); // change to any volumetric mesh above
+    volumeMeshB->computeAttachedSurfaceMesh();
+    auto surfaceMeshB = volumeMeshB->getAttachedSurfaceMesh();
+
+    // Create object B
+    auto objectB = std::make_shared<imstk::VisualObject>("meshObjectVEGA");
+    surfaceMeshB->translate(Vec3d(3, 0, 0));
+    objectB->setVisualGeometry(surfaceMeshB);
+
+    // Add objects to the scene
+    scene->addSceneObject(objectA);
+    scene->addSceneObject(objectB);
+
+    // Run
+    sdk->setCurrentScene("SceneTestMesh");
+    sdk->startSimulation(true);
+
+}
+
 void testVTKTexture()
 {
     // Parse command line arguments
@@ -213,7 +275,7 @@ void testMultiObjectWithTextures()
     auto scene = sdk->createNewScene("multiObjectWithTexturesTest");
 
     // Read surface mesh
-    auto objMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/ETI/resources/OperatingRoom/cloth.obj");
+    auto objMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/ETI/resources/OperatingRoom/cloth.obj");
     auto surfaceMesh = std::dynamic_pointer_cast<imstk::SurfaceMesh>(objMesh);
     surfaceMesh->addTexture(DATA_ROOT_PATH"/ETI/resources/TextureOR/cloth.jpg");
 
@@ -227,7 +289,7 @@ void testMultiObjectWithTextures()
 
     if (secondObject){
         // Read surface mesh1
-        auto objMesh1 = imstk::MeshReader::read(DATA_ROOT_PATH"/ETI/resources/OperatingRoom/bed1.obj");
+        auto objMesh1 = imstk::MeshIO::read(DATA_ROOT_PATH"/ETI/resources/OperatingRoom/bed1.obj");
         auto surfaceMesh1 = std::dynamic_pointer_cast<imstk::SurfaceMesh>(objMesh1);
         if (secondObjectTexture)
             surfaceMesh1->addTexture(DATA_ROOT_PATH"/ETI/resources/TextureOR/bed-1.jpg");
@@ -250,7 +312,7 @@ void testMultiTextures()
     auto scene = sdk->createNewScene("multitexturestest");
 
     // Read surface mesh
-    auto objMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/textures/Fox skull OBJ/fox_skull.obj");
+    auto objMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/textures/Fox skull OBJ/fox_skull.obj");
     auto surfaceMesh = std::dynamic_pointer_cast<imstk::SurfaceMesh>(objMesh);
     surfaceMesh->addTexture(DATA_ROOT_PATH"/textures/Fox skull OBJ/fox_skull_0.jpg",
         "material_0");
@@ -273,8 +335,8 @@ void testMeshCCD()
     auto sdk = std::make_shared<SimulationManager>();
     auto scene = sdk->createNewScene("MeshCCDTest");
 
-    auto mesh1 = imstk::MeshReader::read(DATA_ROOT_PATH"/Spheres/big.vtk");
-    auto mesh2 = imstk::MeshReader::read(DATA_ROOT_PATH"/Spheres/small_0.vtk");
+    auto mesh1 = imstk::MeshIO::read(DATA_ROOT_PATH"/Spheres/big.vtk");
+    auto mesh2 = imstk::MeshIO::read(DATA_ROOT_PATH"/Spheres/small_0.vtk");
 
     // Obj1
     auto obj1 = std::make_shared<CollidingObject>("obj1");
@@ -298,13 +360,13 @@ void testMeshCCD()
     auto t = std::thread([mesh2]
     {
         std::this_thread::sleep_for(std::chrono::seconds(5));
-        auto mesh2_1 = imstk::MeshReader::read(DATA_ROOT_PATH"/Spheres/small_1.vtk");
+        auto mesh2_1 = imstk::MeshIO::read(DATA_ROOT_PATH"/Spheres/small_1.vtk");
         mesh2->setVerticesPositions(mesh2_1->getVertexPositions());
         std::this_thread::sleep_for(std::chrono::seconds(5));
-        auto mesh2_2 = imstk::MeshReader::read(DATA_ROOT_PATH"/Spheres/small_2.vtk");
+        auto mesh2_2 = imstk::MeshIO::read(DATA_ROOT_PATH"/Spheres/small_2.vtk");
         mesh2->setVerticesPositions(mesh2_2->getVertexPositions());
         std::this_thread::sleep_for(std::chrono::seconds(5));
-        auto mesh2_3 = imstk::MeshReader::read(DATA_ROOT_PATH"/Spheres/small_3.vtk");
+        auto mesh2_3 = imstk::MeshIO::read(DATA_ROOT_PATH"/Spheres/small_3.vtk");
         mesh2->setVerticesPositions(mesh2_3->getVertexPositions());
     });
 
@@ -561,7 +623,7 @@ void testCameraController()
     sdk->addModule(client);
 
     // Mesh
-    auto mesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.obj");
+    auto mesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.obj");
     auto meshObject = std::make_shared<imstk::VisualObject>("meshObject");
     meshObject->setVisualGeometry(mesh);
     scene->addSceneObject(meshObject);
@@ -589,15 +651,15 @@ void testReadMesh()
     auto scene = sdk->createNewScene("SceneTestMesh");
 
     // Read surface mesh
-    /*auto objMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.obj");
-    auto plyMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/Cube/models/cube.ply");
-    auto stlMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/Cube/models/cube.stl");
-    auto vtkMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/Cube/models/cube.vtk");
-    auto vtpMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/Cube/models/cube.vtp");*/
+    /*auto objMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.obj");
+    auto plyMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/Cube/models/cube.ply");
+    auto stlMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/Cube/models/cube.stl");
+    auto vtkMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/Cube/models/cube.vtk");
+    auto vtpMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/Cube/models/cube.vtp");*/
 
     // Read volumetricMesh
-    //auto vtkMesh2 = imstk::MeshReader::read(DATA_ROOT_PATH"/AVM/nidus-model/nidus10KTet.vtk");
-    auto vegaMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+    //auto vtkMesh2 = imstk::MeshIO::read(DATA_ROOT_PATH"/AVM/nidus-model/nidus10KTet.vtk");
+    auto vegaMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
 
     // Extract surface mesh
     auto volumeMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(vegaMesh); // change to any volumetric mesh above
@@ -1022,7 +1084,7 @@ void testDeformableBody()
     scene->getCamera()->setPosition(0, 2.0, 15.0);
 
     // b. Load a tetrahedral mesh
-    auto tetMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+    auto tetMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
     if (!tetMesh)
     {
         LOG(WARNING) << "Could not read mesh from file.";
@@ -1062,7 +1124,7 @@ void testDeformableBody()
 
     // Configure dynamic model
     auto dynaModel = std::make_shared<FEMDeformableBodyModel>();
-    dynaModel->configure("asianDragon.config");
+    dynaModel->configure(DATA_ROOT_PATH"/asianDragon/asianDragon.config");
     dynaModel->initialize(volTetMesh);
     auto timeIntegrator = std::make_shared<BackwardEuler>(0.01);// Create and add Backward Euler time integrator
     dynaModel->setTimeIntegrator(timeIntegrator);
@@ -1139,7 +1201,7 @@ void testPbdVolume()
     scene->getCamera()->setPosition(0, 2.0, 15.0);
 
     // b. Load a tetrahedral mesh
-    auto tetMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+    auto tetMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
     if (!tetMesh)
     {
         LOG(WARNING) << "Could not read mesh from file.";
@@ -1300,7 +1362,7 @@ void testPbdCollision()
     scene->getCamera()->setPosition(0, 10.0, 25.0);
 
     // dragon
-    auto tetMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+    auto tetMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
     if (!tetMesh)
     {
         LOG(WARNING) << "Could not read mesh from file.";
@@ -1444,7 +1506,7 @@ void testPbdCollision()
         scene->getCamera()->setPosition(0, 0, 50);
     }
     else if (0){
-        auto tetMesh1 = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+        auto tetMesh1 = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
         if (!tetMesh1)
         {
             LOG(WARNING) << "Could not read mesh from file.";
@@ -1708,9 +1770,9 @@ void testLineMesh()
     else{
         std::string path2obj = DATA_ROOT_PATH"/ETI/resources/Tools/blade2.obj";
 
-        auto collidingMesh = imstk::MeshReader::read(path2obj);
-        auto viusalMesh = imstk::MeshReader::read(path2obj);
-        auto physicsMesh = imstk::MeshReader::read(path2obj);
+        auto collidingMesh = imstk::MeshIO::read(path2obj);
+        auto viusalMesh = imstk::MeshIO::read(path2obj);
+        auto physicsMesh = imstk::MeshIO::read(path2obj);
 
         auto bladeMapP2V = std::make_shared<imstk::OneToOneMap>();
         bladeMapP2V->setMaster(physicsMesh);
@@ -1855,8 +1917,8 @@ void testLineMesh()
         scene->getCamera()->setPosition(0, 0, 50);
     }
     else{
-        //auto tetMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/ETI/resources/Human/tongue.veg");
-        auto tetMesh = imstk::MeshReader::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
+        //auto tetMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/ETI/resources/Human/tongue.veg");
+        auto tetMesh = imstk::MeshIO::read(DATA_ROOT_PATH"/asianDragon/asianDragon.veg");
         if (!tetMesh)
         {
             LOG(WARNING) << "Could not read mesh from file.";