diff --git a/Examples/PBD/PBDClothRemap/CMakeLists.txt b/Examples/PBD/PBDClothRemap/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e86f1ce6919d9a266c8ad6835edeef2eec3160ac
--- /dev/null
+++ b/Examples/PBD/PBDClothRemap/CMakeLists.txt
@@ -0,0 +1,34 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+#  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.
+#
+###########################################################################
+
+project(Example-PBDCloth-Remap)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} pbdClothRemapExample.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/PBD)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} SimulationManager Filtering)
diff --git a/Examples/PBD/PBDClothRemap/pbdClothRemapExample.cpp b/Examples/PBD/PBDClothRemap/pbdClothRemapExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..87882a39cc1b53de57bc692c429fbd7ddc8d3cae
--- /dev/null
+++ b/Examples/PBD/PBDClothRemap/pbdClothRemapExample.cpp
@@ -0,0 +1,244 @@
+/*=========================================================================
+
+   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 "imstkCamera.h"
+#include "imstkKeyboardSceneControl.h"
+#include "imstkLight.h"
+#include "imstkLogger.h"
+#include "imstkMouseSceneControl.h"
+#include "imstkNew.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkRenderMaterial.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkSurfaceMesh.h"
+#include "imstkVisualModel.h"
+#include "imstkVTKViewer.h"
+#include "imstkKeyboardDeviceClient.h"
+#include "imstkSurfaceMeshSubdivide.h"
+
+using namespace imstk;
+
+// Parameters to play with
+const double width  = 10.0;
+const double height = 10.0;
+const int    nRows  = 16;
+const int    nCols  = 16;
+
+///
+/// \brief Creates cloth geometry
+///
+static std::shared_ptr<SurfaceMesh>
+makeClothGeometry(const double width,
+                  const double height,
+                  const int    nRows,
+                  const int    nCols)
+{
+    imstkNew<SurfaceMesh> clothMesh;
+
+    imstkNew<VecDataArray<double, 3>> verticesPtr(nRows * nCols);
+    VecDataArray<double, 3>&          vertices = *verticesPtr.get();
+    const double                      dy       = width / static_cast<double>(nCols - 1);
+    const double                      dx       = height / static_cast<double>(nRows - 1);
+    for (int i = 0; i < nRows; ++i)
+    {
+        for (int j = 0; j < nCols; j++)
+        {
+            vertices[i * nCols + j] = Vec3d(dx * static_cast<double>(i), 1.0, dy * static_cast<double>(j));
+        }
+    }
+
+    // Add connectivity data
+    imstkNew<VecDataArray<int, 3>> indicesPtr;
+    VecDataArray<int, 3>&          indices = *indicesPtr.get();
+    for (int i = 0; i < nRows - 1; ++i)
+    {
+        for (int j = 0; j < nCols - 1; j++)
+        {
+            const int index1 = i * nCols + j;
+            const int index2 = index1 + nCols;
+            const int index3 = index1 + 1;
+            const int index4 = index2 + 1;
+
+            // Interleave [/][\]
+            if (i % 2 ^ j % 2)
+            {
+                indices.push_back(Vec3i(index1, index2, index3));
+                indices.push_back(Vec3i(index4, index3, index2));
+            }
+            else
+            {
+                indices.push_back(Vec3i(index2, index4, index1));
+                indices.push_back(Vec3i(index4, index3, index1));
+            }
+        }
+    }
+
+    clothMesh->initialize(verticesPtr, indicesPtr);
+
+    return clothMesh;
+}
+
+///
+/// \brief Creates cloth object
+///
+static std::shared_ptr<PbdObject>
+makeClothObj(const std::string& name,
+             const double       width,
+             const double       height,
+             const int          nRows,
+             const int          nCols)
+{
+    imstkNew<PbdObject> clothObj(name);
+
+    // Setup the Geometry
+    std::shared_ptr<SurfaceMesh> clothMesh = makeClothGeometry(width, height, nRows, nCols);
+
+    // Setup the Parameters
+    imstkNew<PBDModelConfig> pbdParams;
+    pbdParams->enableConstraint(PbdConstraint::Type::Distance, 1.0e2);
+    pbdParams->enableConstraint(PbdConstraint::Type::Dihedral, 1.0e1);
+    pbdParams->m_fixedNodeIds     = { 0, static_cast<size_t>(nCols) - 1 };
+    pbdParams->m_uniformMassValue = width * height / ((double)nRows * (double)nCols);
+    pbdParams->m_gravity    = Vec3d(0.0, -9.8, 0.0);
+    pbdParams->m_defaultDt  = 0.005;
+    pbdParams->m_iterations = 5;
+
+    // Setup the Model
+    imstkNew<PbdModel> pbdModel;
+    pbdModel->setModelGeometry(clothMesh);
+    pbdModel->configure(pbdParams);
+
+    // Setup the VisualModel
+    imstkNew<RenderMaterial> material;
+    material->setBackFaceCulling(false);
+    material->setDisplayMode(RenderMaterial::DisplayMode::WireframeSurface);
+
+    imstkNew<VisualModel> visualModel(clothMesh);
+    visualModel->setRenderMaterial(material);
+
+    // Setup the Object
+    clothObj->addVisualModel(visualModel);
+    clothObj->setPhysicsGeometry(clothMesh);
+    clothObj->setDynamicalModel(pbdModel);
+
+    return clothObj;
+}
+
+///
+/// \brief This example demonstrates replacement of geometry on a pbd
+/// simualted cloth
+///
+int
+main()
+{
+    // Write log to stdout and file
+    Logger::startLogger();
+
+    // Setup a scene
+    imstkNew<Scene>            scene("PBDCloth");
+    std::shared_ptr<PbdObject> clothObj = nullptr;
+    {
+        clothObj = makeClothObj("Cloth", width, height, nRows, nCols);
+        scene->addSceneObject(clothObj);
+
+        // Light (white)
+        imstkNew<DirectionalLight> whiteLight("whiteLight");
+        whiteLight->setFocalPoint(Vec3d(5.0, -8.0, -5.0));
+        whiteLight->setIntensity(1.0);
+        scene->addLight(whiteLight);
+
+        // Light (red)
+        imstkNew<SpotLight> colorLight("colorLight");
+        colorLight->setPosition(Vec3d(-5.0, -3.0, 5.0));
+        colorLight->setFocalPoint(Vec3d(0.0, -5.0, 5.0));
+        colorLight->setIntensity(100.);
+        colorLight->setColor(Color::Red);
+        colorLight->setSpotAngle(30.0);
+        scene->addLight(colorLight);
+
+        // Adjust camera
+        scene->getActiveCamera()->setFocalPoint(0.0, -5.0, 5.0);
+        scene->getActiveCamera()->setPosition(-15.0, -5.0, 25.0);
+    }
+
+    // Run the simulation
+    {
+        // Setup a viewer to render in its own thread
+        imstkNew<VTKViewer> viewer("Viewer");
+        viewer->setActiveScene(scene);
+
+        // Setup a scene manager to advance the scene in its own thread
+        imstkNew<SceneManager> sceneManager("Scene Manager");
+        sceneManager->setActiveScene(scene);
+        viewer->addChildThread(sceneManager); // SceneManager will start/stop with viewer
+
+        // Add mouse and keyboard controls to the viewer
+        {
+            imstkNew<MouseSceneControl> mouseControl(viewer->getMouseDevice());
+            mouseControl->setSceneManager(sceneManager);
+            viewer->addControl(mouseControl);
+
+            imstkNew<KeyboardSceneControl> keyControl(viewer->getKeyboardDevice());
+            keyControl->setSceneManager(sceneManager);
+            keyControl->setViewer(viewer);
+            viewer->addControl(keyControl);
+        }
+
+        // Queue keypress to be called after scene thread
+        queueConnect<KeyPressEvent>(viewer->getKeyboardDevice(), EventType::KeyEvent, sceneManager, [&](KeyPressEvent* e)
+        {
+            // When i is pressed replace the PBD cloth with a subdivided one
+            if (e->m_key == 'i' && e->m_keyPressType == KEY_PRESS)
+            {
+                // This has a number of issues that make it not physically realistic
+                // - Mass is not conservative when interpolated from subdivide
+                // - Constraint resting lengths are not correctly reinited
+                std::shared_ptr<SurfaceMesh> clothMesh = std::dynamic_pointer_cast<SurfaceMesh>(clothObj->getPhysicsGeometry());
+                imstkNew<SurfaceMeshSubdivide> subdiv;
+                subdiv->setInputMesh(clothMesh);
+                subdiv->setNumberOfSubdivisions(1);
+                subdiv->setSubdivisionType(SurfaceMeshSubdivide::Type::LINEAR);
+                subdiv->update();
+                std::shared_ptr<SurfaceMesh> newClothMesh = subdiv->getOutputMesh();
+
+                // RenderDelegates cannot visually have entire geometries swapped yet, so even though we could just set the geometry
+                // on the model, you would not visually see it. Instead we replace the vertex and index buffers of the existing one
+                // Another issue here is that initial geometry is not remapped so reset will not reset to undeformed config
+                clothMesh->setInitialVertexPositions(std::make_shared<VecDataArray<double, 3>>(*newClothMesh->getVertexPositions()));
+                clothMesh->setVertexPositions(newClothMesh->getVertexPositions());
+                clothMesh->setTriangleIndices(newClothMesh->getTriangleIndices());
+                clothMesh->setVertexAttribute("Velocities", newClothMesh->getVertexAttribute("Velocities"));
+                clothMesh->modified();
+
+                // Re-setup the constraints on the object
+                clothObj->initialize();
+            }
+            });
+
+        // Start viewer running, scene as paused
+        sceneManager->requestStatus(ThreadStatus::Paused);
+        viewer->start();
+    }
+
+    return 0;
+}
diff --git a/Examples/Rendering/RenderingExample.cpp b/Examples/Rendering/RenderingExample.cpp
index 4d3b7d2fc65e65783050173d5d99e0213e455a39..6af09bf96c3483930aa50b2563ecf4ba46766e18 100644
--- a/Examples/Rendering/RenderingExample.cpp
+++ b/Examples/Rendering/RenderingExample.cpp
@@ -76,6 +76,7 @@ main()
         material->addTexture(headDiffuseTexture);
         material->addTexture(headNormalTexture);
         material->addTexture(headAoTexture);
+        material->setRecomputeVertexNormals(false);
 
         imstkNew<VisualModel> surfMeshModel(surfaceMesh);
         surfMeshModel->setRenderMaterial(material);
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index ce0cd14c46f5c54398d5f50777ae9edf0f4166be..3b2daf195abf547c7d42f2ac69b9b804c22edbfa 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -42,8 +42,6 @@ PbdModel::PbdModel() : DynamicalModel(DynamicalModelType::PositionBasedDynamics)
     m_mass(std::make_shared<DataArray<double>>()),
     m_invMass(std::make_shared<DataArray<double>>()),
     m_fixedNodeInvMass(std::make_shared<std::unordered_map<size_t, double>>()),
-    m_constraints(std::make_shared<PBDConstraintVector>()),
-    m_partitionedConstraints(std::make_shared<std::vector<PBDConstraintVector>>()),
     m_parameters(std::make_shared<PBDModelConfig>())
 {
     m_validGeometryTypes = {
@@ -103,87 +101,91 @@ bool
 PbdModel::initialize()
 {
     LOG_IF(FATAL, (!this->getModelGeometry())) << "Model geometry is not yet set! Cannot initialize without model geometry.";
+    bool bOK = true; // Return immediately if some constraint failed to initialize
 
     initState();
 
-    bool bOK = true; // Return immediately if some constraint failed to initialize
-
-    // Setup the default pbd solver if none exists
-    if (m_pbdSolver == nullptr)
+    // Initialize constraints
     {
-        m_pbdSolver = std::make_shared<PbdSolver>();
-        m_pbdSolver->setIterations(m_parameters->m_iterations);
-        m_pbdSolver->setSolverType(m_parameters->m_solverType);
-    }
-    m_pbdSolver->setPositions(getCurrentState()->getPositions());
-    m_pbdSolver->setInvMasses(getInvMasses());
-    m_pbdSolver->setConstraints(getConstraints());
-    m_pbdSolver->setPartitionedConstraints(getPartitionedConstraints());
-    m_pbdSolver->setTimeStep(m_parameters->m_dt);
+        m_constraints = std::make_shared<PBDConstraintVector>();
 
-    // Initialize FEM constraints
-    for (auto& constraint: m_parameters->m_FEMConstraints)
-    {
-        computeElasticConstants();
-        if (!initializeFEMConstraints(constraint.second))
+        // Initialize FEM constraints
+        for (auto& constraint : m_parameters->m_FEMConstraints)
         {
-            return false;
+            computeElasticConstants();
+            if (!initializeFEMConstraints(constraint.second))
+            {
+                return false;
+            }
         }
-    }
 
-    // Initialize other constraints
-    for (auto& constraint: m_parameters->m_regularConstraints)
-    {
-        if (m_parameters->m_solverType == PbdConstraint::SolverType::PBD && constraint.second > 1.0)
-        {
-            LOG(WARNING) << "for PBD, k should be between [0, 1]";
-        }
-        else if (m_parameters->m_solverType == PbdConstraint::SolverType::xPBD && constraint.second <= 1.0)
+        // Initialize other constraints
+        for (auto& constraint : m_parameters->m_regularConstraints)
         {
-            LOG(WARNING) << "for xPBD, k is Young's Modulu, and should be much larger than 1";
-        }
+            if (m_parameters->m_solverType == PbdConstraint::SolverType::PBD && constraint.second > 1.0)
+            {
+                LOG(WARNING) << "for PBD, k should be between [0, 1]";
+            }
+            else if (m_parameters->m_solverType == PbdConstraint::SolverType::xPBD && constraint.second <= 1.0)
+            {
+                LOG(WARNING) << "for xPBD, k is Young's Modulu, and should be much larger than 1";
+            }
 
-        if (!bOK)
-        {
-            return false;
-        }
-        switch (constraint.first)
-        {
-        case PbdConstraint::Type::Volume:
-            bOK = initializeVolumeConstraints(constraint.second);
-            break;
+            if (!bOK)
+            {
+                return false;
+            }
+            switch (constraint.first)
+            {
+            case PbdConstraint::Type::Volume:
+                bOK = initializeVolumeConstraints(constraint.second);
+                break;
 
-        case PbdConstraint::Type::Distance:
-            bOK = initializeDistanceConstraints(constraint.second);
-            break;
+            case PbdConstraint::Type::Distance:
+                bOK = initializeDistanceConstraints(constraint.second);
+                break;
 
-        case PbdConstraint::Type::Area:
-            bOK = initializeAreaConstraints(constraint.second);
-            break;
+            case PbdConstraint::Type::Area:
+                bOK = initializeAreaConstraints(constraint.second);
+                break;
 
-        case PbdConstraint::Type::Bend:
-            bOK = initializeBendConstraints(constraint.second);
-            break;
+            case PbdConstraint::Type::Bend:
+                bOK = initializeBendConstraints(constraint.second);
+                break;
 
-        case PbdConstraint::Type::Dihedral:
-            bOK = initializeDihedralConstraints(constraint.second);
-            break;
+            case PbdConstraint::Type::Dihedral:
+                bOK = initializeDihedralConstraints(constraint.second);
+                break;
 
-        case PbdConstraint::Type::ConstantDensity:
-            bOK = initializeConstantDensityConstraint(constraint.second);
-            break;
+            case PbdConstraint::Type::ConstantDensity:
+                bOK = initializeConstantDensityConstraint(constraint.second);
+                break;
+
+            default:
+                LOG(FATAL) << "Invalid constraint type";
+            }
+        }
 
-        default:
-            LOG(FATAL) << "Invalid constraint type";
+        // Partition constraints for parallel computation
+        if (!m_partitioned)
+        {
+            this->partitionConstraints();
+            m_partitioned = true;
         }
     }
 
-    // Partition constraints for parallel computation
-    if (!m_partitioned)
+    // Setup the default pbd solver if none exists
+    if (m_pbdSolver == nullptr)
     {
-        this->partitionConstraints();
-        m_partitioned = true;
+        m_pbdSolver = std::make_shared<PbdSolver>();
+        m_pbdSolver->setIterations(m_parameters->m_iterations);
+        m_pbdSolver->setSolverType(m_parameters->m_solverType);
     }
+    m_pbdSolver->setPositions(getCurrentState()->getPositions());
+    m_pbdSolver->setInvMasses(getInvMasses());
+    m_pbdSolver->setConstraints(getConstraints());
+    m_pbdSolver->setPartitionedConstraints(getPartitionedConstraints());
+    m_pbdSolver->setTimeStep(m_parameters->m_dt);
 
     this->setTimeStepSizeType(m_timeStepSizeType);
 
@@ -235,8 +237,23 @@ PbdModel::initState()
         }
     }
 
-    // Define velocities and accelerations on the geometry
-    m_mesh->setVertexAttribute("Velocities", m_currentState->getVelocities());
+    // Initialize Velocities
+    {
+        // If the input mesh has per vertex velocities, use those
+        std::shared_ptr<AbstractDataArray> velocities = m_mesh->getVertexAttribute("Velocities");
+        if (velocities != nullptr && velocities->getNumberOfComponents() == 3 && velocities->getScalarType() == IMSTK_DOUBLE
+            && std::dynamic_pointer_cast<VecDataArray<double, 3>>(velocities)->size() == numParticles)
+        {
+            m_currentState->setVelocities(std::dynamic_pointer_cast<VecDataArray<double, 3>>(velocities));
+        }
+        // If not, put existing (0 initialized velocities) on mesh
+        else
+        {
+            m_mesh->setVertexAttribute("Velocities", m_currentState->getVelocities());
+        }
+    }
+
+    // Define accelerations on the geometry
     m_mesh->setVertexAttribute("Accelerations", m_currentState->getAccelerations());
 
     // Overwrite some masses for specified fixed points
@@ -560,6 +577,8 @@ PbdModel::initializeConstantDensityConstraint(const double stiffness)
 void
 PbdModel::partitionConstraints(const bool print)
 {
+    m_partitionedConstraints = std::make_shared<std::vector<PBDConstraintVector>>();
+
     // Form the map { vertex : list_of_constraints_involve_vertex }
     PBDConstraintVector& allConstraints = *m_constraints;
 
diff --git a/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp b/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
index d5e86be8afe2b3052cc6dc3c5060f892651b827d..089b0107676f4300c55cfa367704248364b076a7 100644
--- a/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
+++ b/Source/DynamicalModels/ObjectStates/imstkPbdState.cpp
@@ -29,5 +29,7 @@ PbdState::setState(std::shared_ptr<PbdState> rhs)
     *m_pos = *rhs->getPositions();
     *m_vel = *rhs->getVelocities();
     *m_acc = *rhs->getAccelerations();
+
+    m_pos->modified();
 }
 } // imstk
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp b/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
index f76a23c8984dea960cba669a2bb91ab43ffa127b..a492b7ed8ae870d307128854b7c56d1fdd6e70a7 100644
--- a/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
+++ b/Source/DynamicalModels/ObjectStates/imstkSPHState.cpp
@@ -65,6 +65,8 @@ SPHState::setState(std::shared_ptr<SPHState> rhs)
     m_NeighborLists   = rhs->getFluidNeighborLists();
     m_BDNeighborLists = rhs->getBoundaryNeighborLists();
     m_NeighborInfo    = rhs->getNeighborInfo();
+
+    m_positions->modified();
 }
 
 size_t
diff --git a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKHexahedralMeshRenderDelegate.cpp b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKHexahedralMeshRenderDelegate.cpp
index e12a8f1c8ea221e520dfc818e9c6b810728d1e6d..df0e75b5b332f99c885f3d319f8a8b7d6fbd153f 100644
--- a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKHexahedralMeshRenderDelegate.cpp
+++ b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKHexahedralMeshRenderDelegate.cpp
@@ -162,7 +162,7 @@ VTKHexahedralMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
         m_indices = geometry->getHexahedraIndices();
         {
             // Copy cells
-            m_cellArray = vtkSmartPointer<vtkCellArray>::New();
+            m_cellArray->Reset();
             vtkIdType cell[8];
             for (const auto& t : *m_indices)
             {
@@ -172,8 +172,7 @@ VTKHexahedralMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
                 }
                 m_cellArray->InsertNextCell(8, cell);
             }
-            m_mesh->SetCells(VTK_HEXAHEDRON, m_cellArray);
-            m_mesh->Modified();
+            m_cellArray->Modified();
         }
         vertexOrIndexBufferChanged = true;
     }
diff --git a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKLineMeshRenderDelegate.cpp b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKLineMeshRenderDelegate.cpp
index e8caf76e4b2fc45a02bf778235ac93ccfa7c37ff..70a058f8035c7903b8fa00498f00abc522afe1db 100644
--- a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKLineMeshRenderDelegate.cpp
+++ b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKLineMeshRenderDelegate.cpp
@@ -181,7 +181,7 @@ VTKLineMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
         m_indices = geometry->getLinesIndices();
         {
             // Copy cells
-            m_cellArray = vtkSmartPointer<vtkCellArray>::New();
+            m_cellArray->Reset();
             vtkIdType cell[2];
             for (const auto& t : *m_indices)
             {
@@ -191,8 +191,7 @@ VTKLineMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
                 }
                 m_cellArray->InsertNextCell(2, cell);
             }
-            m_polydata->SetPolys(m_cellArray);
-            m_polydata->Modified();
+            m_cellArray->Modified();
         }
     }
 }
diff --git a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.cpp b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.cpp
index 1678f5199acb512221bf249a8d2cf525c05e6882..328d8b0376838882a08d83cba9101146689d6fcd 100644
--- a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.cpp
+++ b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.cpp
@@ -44,19 +44,19 @@ VTKSurfaceMeshRenderDelegate::VTKSurfaceMeshRenderDelegate(std::shared_ptr<Visua
     m_mappedVertexArray(vtkSmartPointer<vtkDoubleArray>::New()),
     m_mappedNormalArray(vtkSmartPointer<vtkDoubleArray>::New())
 {
-    auto geometry = std::static_pointer_cast<SurfaceMesh>(m_visualModel->getGeometry());
-    geometry->computeVertexNeighborTriangles();
+    m_geometry = std::static_pointer_cast<SurfaceMesh>(m_visualModel->getGeometry());
+    m_geometry->computeVertexNeighborTriangles();
 
     // Get our own handles to these in case the geometry changes them
-    m_vertices = geometry->getVertexPositions();
-    m_indices  = geometry->getTriangleIndices();
+    m_vertices = m_geometry->getVertexPositions();
+    m_indices  = m_geometry->getTriangleIndices();
 
     // Map vertices to VTK point data
     if (m_vertices != nullptr)
     {
         m_mappedVertexArray = vtkDoubleArray::SafeDownCast(GeometryUtils::coupleVtkDataArray(m_vertices));
         auto points = vtkSmartPointer<vtkPoints>::New();
-        points->SetNumberOfPoints(geometry->getNumVertices());
+        points->SetNumberOfPoints(m_geometry->getNumVertices());
         points->SetData(m_mappedVertexArray);
         m_polydata->SetPoints(points);
     }
@@ -78,32 +78,32 @@ VTKSurfaceMeshRenderDelegate::VTKSurfaceMeshRenderDelegate(std::shared_ptr<Visua
     }
 
     // Map vertex scalars if it has them
-    if (geometry->getVertexScalars() != nullptr)
+    if (m_geometry->getVertexScalars() != nullptr)
     {
-        m_mappedVertexScalarArray = GeometryUtils::coupleVtkDataArray(geometry->getVertexScalars());
+        m_mappedVertexScalarArray = GeometryUtils::coupleVtkDataArray(m_geometry->getVertexScalars());
         m_polydata->GetPointData()->SetScalars(m_mappedVertexScalarArray);
     }
 
     // Map cell scalars if it has them
-    if (geometry->getCellScalars() != nullptr)
+    if (m_geometry->getCellScalars() != nullptr)
     {
-        m_mappedCellScalarArray = GeometryUtils::coupleVtkDataArray(geometry->getCellScalars());
+        m_mappedCellScalarArray = GeometryUtils::coupleVtkDataArray(m_geometry->getCellScalars());
         m_polydata->GetCellData()->SetScalars(m_mappedCellScalarArray);
     }
 
     // Map normals, if none provided compute per vertex normals
-    if (geometry->getVertexNormals() == nullptr)
+    if (m_geometry->getVertexNormals() == nullptr)
     {
-        geometry->computeVertexNormals();
+        m_geometry->computeVertexNormals();
     }
-    m_mappedNormalArray = vtkDoubleArray::SafeDownCast(GeometryUtils::coupleVtkDataArray(geometry->getVertexNormals()));
+    m_mappedNormalArray = vtkDoubleArray::SafeDownCast(GeometryUtils::coupleVtkDataArray(m_geometry->getVertexNormals()));
     m_polydata->GetPointData()->SetNormals(m_mappedNormalArray);
 
     // Map TCoords
-    if (geometry->getVertexTCoords() != nullptr)
+    if (m_geometry->getVertexTCoords() != nullptr)
     {
-        m_mappedTCoordsArray = vtkFloatArray::SafeDownCast(GeometryUtils::coupleVtkDataArray(geometry->getVertexTCoords()));
-        m_mappedTCoordsArray->SetName(geometry->getActiveVertexTCoords().c_str());
+        m_mappedTCoordsArray = vtkFloatArray::SafeDownCast(GeometryUtils::coupleVtkDataArray(m_geometry->getVertexTCoords()));
+        m_mappedTCoordsArray->SetName(m_geometry->getActiveVertexTCoords().c_str());
         m_polydata->GetPointData()->SetTCoords(m_mappedTCoordsArray);
 
         // Map Tangents
@@ -117,10 +117,12 @@ VTKSurfaceMeshRenderDelegate::VTKSurfaceMeshRenderDelegate(std::shared_ptr<Visua
     }
 
     // When geometry is modified, update data source, mostly for when an entirely new array/buffer was set
-    queueConnect<Event>(geometry, EventType::Modified, this, &VTKSurfaceMeshRenderDelegate::geometryModified);
+    queueConnect<Event>(m_geometry, EventType::Modified, this, &VTKSurfaceMeshRenderDelegate::geometryModified);
 
     // When the vertex buffer internals are modified, ie: a single or N elements
-    queueConnect<Event>(geometry->getVertexPositions(), EventType::Modified, this, &VTKSurfaceMeshRenderDelegate::vertexDataModified);
+    queueConnect<Event>(m_geometry->getVertexPositions(), EventType::Modified, this, &VTKSurfaceMeshRenderDelegate::vertexDataModified);
+
+    // When the index buffer internals are modified,
 
     // Setup mapper
     {
@@ -231,6 +233,7 @@ VTKSurfaceMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
             m_mappedVertexArray->SetNumberOfComponents(3);
             m_mappedVertexArray->SetArray(reinterpret_cast<double*>(m_vertices->getPointer()), m_vertices->size() * 3, 1);
         }
+        m_polydata->GetPoints()->SetNumberOfPoints(m_vertices->size());
         //vertexOrIndexBufferChanged = true;
     }
 
@@ -243,7 +246,7 @@ VTKSurfaceMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
         m_indices = geometry->getTriangleIndices();
         {
             // Copy cells
-            m_cellArray = vtkSmartPointer<vtkCellArray>::New();
+            m_cellArray->Reset();
             vtkIdType cell[3];
             for (const auto& t : *m_indices)
             {
@@ -253,8 +256,7 @@ VTKSurfaceMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
                 }
                 m_cellArray->InsertNextCell(3, cell);
             }
-            m_polydata->SetPolys(m_cellArray);
-            m_polydata->Modified();
+            m_cellArray->Modified();
         }
         //vertexOrIndexBufferChanged = true;
     }
diff --git a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.h b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.h
index 96d3fa8b8ccd2b7f296dfb108d780cf0e66857b2..d8e212ee644e127b6c80328bee626ae3cfea4e44 100644
--- a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.h
+++ b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKSurfaceMeshRenderDelegate.h
@@ -77,6 +77,7 @@ protected:
     void geometryModified(Event* e);
 
 protected:
+    std::shared_ptr<SurfaceMesh> m_geometry;
     std::shared_ptr<VecDataArray<double, 3>> m_vertices;
     std::shared_ptr<VecDataArray<double, 3>> m_normals;
     std::shared_ptr<VecDataArray<int, 3>>    m_indices;
diff --git a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKTetrahedralMeshRenderDelegate.cpp b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKTetrahedralMeshRenderDelegate.cpp
index cc78f31f60e59afeefad31c7d6db9a4c7f77b010..accf96d77ecdc8a6ee0b39acda064910b4f74962 100644
--- a/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKTetrahedralMeshRenderDelegate.cpp
+++ b/Source/Rendering/VTKRenderer/RenderDelegate/imstkVTKTetrahedralMeshRenderDelegate.cpp
@@ -181,7 +181,7 @@ VTKTetrahedralMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
         m_indices = geometry->getTetrahedraIndices();
         {
             // Copy cells
-            m_cellArray = vtkSmartPointer<vtkCellArray>::New();
+            m_cellArray->Reset();
             vtkIdType cell[3];
             for (const auto& t : *m_indices)
             {
@@ -191,8 +191,7 @@ VTKTetrahedralMeshRenderDelegate::geometryModified(Event* imstkNotUsed(e))
                 }
                 m_cellArray->InsertNextCell(3, cell);
             }
-            m_mesh->SetCells(VTK_TETRA, m_cellArray);
-            m_mesh->Modified();
+            m_cellArray->Modified();
         }
     }
 }