diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index aabe9cfe8b57c0192b2cf8e5009aca35c0d6a405..6ff7dd40ae4d2485932423f8727f10c6756af838 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -264,21 +264,9 @@ PbdModel::initializeDistanceConstraints(const double stiffness)
     }
     else if (m_mesh->getType() == Geometry::Type::SurfaceMesh)
     {
-        // To remove duplicate edges convert to a line mesh (indices are then messed up)
-        std::shared_ptr<LineMesh>      lineMesh = GeometryUtils::SurfaceMeshToLineMesh(std::static_pointer_cast<SurfaceMesh>(m_mesh));
-        MeshIO::write(lineMesh, "C:/Users/andrew.wilson/Desktop/test.vtk");
-        const auto&                    elements = lineMesh->getLinesVertices();
-        const auto                     nV       = lineMesh->getNumVertices();
-        std::vector<std::vector<bool>> E(nV, std::vector<bool>(nV, 1));
-
-        for (size_t k = 0; k < elements.size(); ++k)
-        {
-            auto& seg = elements[k];
-            addConstraint(E, seg[0], seg[1]);
-        }
-        /*const auto& triMesh  = std::static_pointer_cast<SurfaceMesh>(m_mesh);
-        const auto& elements = triMesh->getTrianglesVertices();
-        const auto        nV = triMesh->getNumVertices();
+        const auto&                    triMesh  = std::static_pointer_cast<SurfaceMesh>(m_mesh);
+        const auto&                    elements = triMesh->getTrianglesVertices();
+        const auto                     nV       = triMesh->getNumVertices();
         std::vector<std::vector<bool>> E(nV, std::vector<bool>(nV, 1));
 
         for (size_t k = 0; k < elements.size(); ++k)
@@ -287,8 +275,7 @@ PbdModel::initializeDistanceConstraints(const double stiffness)
             addConstraint(E, tri[0], tri[1]);
             addConstraint(E, tri[0], tri[2]);
             addConstraint(E, tri[1], tri[2]);
-        }*/
-
+        }
     }
     else if (m_mesh->getType() == Geometry::Type::LineMesh)
     {
diff --git a/Source/Geometry/Reader/imstkVTKMeshIO.cpp b/Source/Geometry/Reader/imstkVTKMeshIO.cpp
index c985888c50d7a16af5bf54fbed5d4643720d4ee9..4906ec92901a347be69a0f14b583686dba0f5a89 100644
--- a/Source/Geometry/Reader/imstkVTKMeshIO.cpp
+++ b/Source/Geometry/Reader/imstkVTKMeshIO.cpp
@@ -35,8 +35,7 @@
 #include "vtkSTLWriter.h"
 #include "vtkFloatArray.h"
 #include "vtkTriangleFilter.h"
-#include "vtkDICOMImageReader.h"
-#include "vtkNrrdReader.h"
+#include "vtkPolyDataWriter.h"
 
 #include "g3log/g3log.hpp"
 
@@ -111,6 +110,8 @@ VTKMeshIO::write(const std::shared_ptr<PointSet> imstkMesh, const std::string& f
             return VTKMeshIO::writeVtkPolyData<vtkSTLWriter>(sMesh, filePath);
         case MeshFileType::PLY:
             return VTKMeshIO::writeVtkPolyData<vtkPLYWriter>(sMesh, filePath);
+        case MeshFileType::VTK:
+            return VTKMeshIO::writeVtkPolyData<vtkPolyDataWriter>(sMesh, filePath);
         default:
             LOG(WARNING) << "VTKMeshIO::write error: file type not supported for surface mesh.";
             return false;
@@ -121,6 +122,8 @@ VTKMeshIO::write(const std::shared_ptr<PointSet> imstkMesh, const std::string& f
         switch (meshType)
         {
         case MeshFileType::VTK:
+            return VTKMeshIO::writeVtkPolyData<vtkPolyDataWriter>(lMesh, filePath);
+        case MeshFileType::VTP:
             return VTKMeshIO::writeVtkPolyData<vtkXMLPolyDataWriter>(lMesh, filePath);
         default:
             LOG(WARNING) << "vtkMeshIO::write error: file type not supported for line mesh.";
@@ -196,6 +199,7 @@ bool
 VTKMeshIO::writeVtkPolyData(const std::shared_ptr<LineMesh> imstkMesh, const std::string& filePath)
 {
     vtkSmartPointer<vtkPolyData> vtkMesh = GeometryUtils::convertLineMeshToVtkPolyData(imstkMesh);
+    vtkMesh->PrintSelf(std::cout, vtkIndent());
     if (!vtkMesh)  //conversion unsuccessful
     {
         return false;
diff --git a/Source/Geometry/imstkGeometryUtilities.cpp b/Source/Geometry/imstkGeometryUtilities.cpp
index 2f9819e269cac7026b70d6151c0bcd51afaceb6a..719b74614afc2f7913c2f548017e6edfa0a4a59d 100644
--- a/Source/Geometry/imstkGeometryUtilities.cpp
+++ b/Source/Geometry/imstkGeometryUtilities.cpp
@@ -1,4 +1,15 @@
 #include "imstkGeometryUtilities.h"
+#include "imstkHexahedralMesh.h"
+#include "imstkLineMesh.h"
+#include "imstkTetrahedralMesh.h"
+#include <vtkAppendPolyData.h>
+#include <vtkExtractEdges.h>
+#include <vtkLinearSubdivisionFilter.h>
+#include <vtkLoopSubdivisionFilter.h>
+#include <vtkPointData.h>
+#include <vtkSmoothPolyDataFilter.h>
+#include <vtkTriangleFilter.h>
+#include <vtkUnstructuredGrid.h>
 
 namespace imstk
 {
@@ -7,7 +18,7 @@ GeometryUtils::convertVtkPolyDataToSurfaceMesh(vtkSmartPointer<vtkPolyData> vtkM
 {
     if (!vtkMesh)
     {
-        LOG(FATAL) << "VTKMeshIO::convertVtkPolyDataToSurfaceMesh error: could not read with VTK reader.";
+        LOG(FATAL) << "GeometryUtils::convertVtkPolyDataToSurfaceMesh error: could not convert vtkPolyData to imstkSurfaceMesh.";
         return nullptr;
     }
 
@@ -45,7 +56,7 @@ GeometryUtils::convertVtkPolyDataToLineMesh(vtkSmartPointer<vtkPolyData> vtkMesh
 {
     if (!vtkMesh)
     {
-        LOG(FATAL) << "VTKMeshIO::convertVtkPolyDataToSurfaceMesh error: could not read with VTK reader.";
+        LOG(FATAL) << "GeometryUtils::convertVtkPolyDataToSurfaceMesh error: could not convert vtkPolyData to imstkLineMesh.";
         return nullptr;
     }
 
@@ -79,7 +90,7 @@ GeometryUtils::convertSurfaceMeshToVtkPolyData(std::shared_ptr<SurfaceMesh> imst
     vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
     copyCellsToVtk<3>(imstkMesh->getTrianglesVertices(), polys.Get());
 
-    vtkPolyData* polydata = vtkPolyData::New();
+    vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
     polydata->SetPoints(points);
     polydata->SetPolys(polys);
     return polydata;
@@ -94,7 +105,7 @@ GeometryUtils::convertLineMeshToVtkPolyData(std::shared_ptr<LineMesh> imstkMesh)
     vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
     copyCellsToVtk<2>(imstkMesh->getLinesVertices(), polys.Get());
 
-    vtkPolyData* polydata = vtkPolyData::New();
+    vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
     polydata->SetPoints(points);
     polydata->SetPolys(polys);
     return polydata;
@@ -169,6 +180,7 @@ GeometryUtils::convertVtkUnstructuredGridToVolumetricMesh(vtkUnstructuredGrid* v
     }
 }
 
+
 void
 GeometryUtils::copyVerticesFromVtk(vtkPoints* points, StdVectorOfVec3d& vertices)
 {
@@ -235,7 +247,7 @@ GeometryUtils::copyCellsFromVtk(vtkCellArray* vtkCells, std::vector<std::array<s
 
     cells.reserve(vtkCells->GetNumberOfCells());
     vtkCells->InitTraversal();
-    auto                    vtkCell = vtkSmartPointer<vtkIdList>::New();
+    vtkNew<vtkIdList> vtkCell;
     std::array<size_t, dim> cell;
     while (vtkCells->GetNextCell(vtkCell))
     {
@@ -281,10 +293,11 @@ GeometryUtils::copyPointDataFromVtk(vtkPointData* pointData, std::map<std::strin
     }
 }
 
+
 std::shared_ptr<SurfaceMesh>
 GeometryUtils::appendSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh1, std::shared_ptr<SurfaceMesh> surfaceMesh2)
 {
-    vtkSmartPointer<vtkAppendPolyData> filter = vtkSmartPointer<vtkAppendPolyData>::New();
+    vtkNew<vtkAppendPolyData> filter;
     filter->AddInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh1));
     filter->AddInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh2));
     filter->Update();
@@ -292,11 +305,54 @@ GeometryUtils::appendSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh1, std:
 }
 
 std::shared_ptr<LineMesh>
-GeometryUtils::SurfaceMeshToLineMesh(std::shared_ptr<SurfaceMesh> surfMesh)
+GeometryUtils::surfaceMeshToLineMesh(std::shared_ptr<SurfaceMesh> surfaceMesh)
+{
+    vtkNew<vtkExtractEdges> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->Update();
+    vtkNew<vtkTriangleFilter> triangleFilter;
+    triangleFilter->SetInputData(filter->GetOutput());
+    triangleFilter->Update();
+    return convertVtkPolyDataToLineMesh(triangleFilter->GetOutput());
+}
+
+std::shared_ptr<SurfaceMesh>
+GeometryUtils::smoothSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh,
+    int numberOfIterations, double relaxationFactor,
+    double convergence, double featureAngle,
+    double edgeAngle, bool featureEdgeSmoothing,
+    bool boundarySmoothing)
+{
+    vtkNew<vtkSmoothPolyDataFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfIterations(numberOfIterations);
+    filter->SetRelaxationFactor(relaxationFactor);
+    filter->SetConvergence(convergence);
+    filter->SetFeatureAngle(featureAngle);
+    filter->SetEdgeAngle(edgeAngle);
+    filter->SetFeatureEdgeSmoothing(featureEdgeSmoothing);
+    filter->SetBoundarySmoothing(boundarySmoothing);
+    filter->Update();
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+
+std::shared_ptr<SurfaceMesh>
+GeometryUtils::linearSubdivideSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, int numberOfSubdivisions)
+{
+    vtkNew<vtkLinearSubdivisionFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfSubdivisions(numberOfSubdivisions);
+    filter->Update();
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
+}
+
+std::shared_ptr<SurfaceMesh>
+GeometryUtils::loopSubdivideSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, int numberOfSubdivisions)
 {
-    vtkSmartPointer<vtkExtractEdges> filter = vtkSmartPointer<vtkExtractEdges>::New();
-    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfMesh));
+    vtkNew<vtkLoopSubdivisionFilter> filter;
+    filter->SetInputData(convertSurfaceMeshToVtkPolyData(surfaceMesh));
+    filter->SetNumberOfSubdivisions(numberOfSubdivisions);
     filter->Update();
-    return convertVtkPolyDataToLineMesh(filter->GetOutput());
+    return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
 }
 }
\ No newline at end of file
diff --git a/Source/Geometry/imstkGeometryUtilities.h b/Source/Geometry/imstkGeometryUtilities.h
index 9d6c7d230da49f07b5e01920528101ad26e501a6..e9597f36f3be8b84dfd1693f57d68a38115e4f09 100644
--- a/Source/Geometry/imstkGeometryUtilities.h
+++ b/Source/Geometry/imstkGeometryUtilities.h
@@ -21,20 +21,23 @@
 
 #pragma once
 
-#include "imstkLineMesh.h"
-#include "imstkSurfaceMesh.h"
-#include "imstkVTKMeshIO.h"
+#include "imstkMath.h"
 #include <memory>
-#include <vtkAppendPolyData.h>
-#include <vtkCellArray.h>
-#include <vtkCleanPolyData.h>
-#include <vtkPoints.h>
-#include <vtkPolyData.h>
 #include <vtkSmartPointer.h>
-#include <vtkExtractEdges.h>
+
+class vtkCellArray;
+class vtkPolyData;
+class vtkPointData;
+class vtkPoints;
+class vtkUnstructuredGrid;
 
 namespace imstk
 {
+class HexahedralMesh;
+class LineMesh;
+class SurfaceMesh;
+class TetrahedralMesh;
+class VolumetricMesh;
 
 class GeometryUtils
 {
@@ -108,8 +111,27 @@ public:
     static std::shared_ptr<SurfaceMesh> appendSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh1, std::shared_ptr<SurfaceMesh> surfaceMesh2);
 
     ///
-    /// \brief Converts an imstk SurfaceMesh to a LineMesh, removing duplicate edges
+    /// \brief Converts an imstk SurfaceMesh to a LineMesh, removing duplicate edges. Cell indices not preserved
+    ///
+    static std::shared_ptr<LineMesh> surfaceMeshToLineMesh(std::shared_ptr<SurfaceMesh> surfaceMesh);
+
+    ///
+    /// \brief Smooths a SurfaceMesh using laplacian smoothening
+    ///
+    static std::shared_ptr<SurfaceMesh> smoothSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh,
+        int numberOfIterations = 20, double relaxationFactor = 0.01,
+        double convergence = 0.0, double featureAngle = 45.0,
+        double edgeAngle = 15.0, bool featureEdgeSmoothing = false,
+        bool boundarySmoothing = true);
+
+    ///
+    /// \brief Subidivdes a SurfaceMesh using linear subdivision
+    ///
+    static std::shared_ptr<SurfaceMesh> linearSubdivideSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, int numberOfSubdivisions);
+
+    ///
+    /// \brief Subidivides a SurfaceMesh using loop subdivision
     ///
-    static std::shared_ptr<LineMesh> SurfaceMeshToLineMesh(std::shared_ptr<SurfaceMesh> surfMesh);
+    static std::shared_ptr<SurfaceMesh> loopSubdivideSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, int numberOfSubdivisions);
 };
 }