diff --git a/Examples/PBDCutting/CMakeLists.txt b/Examples/PBDCutting/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a2aff794533c4c81eed61e23e58e45c4b97ec27f
--- /dev/null
+++ b/Examples/PBDCutting/CMakeLists.txt
@@ -0,0 +1,37 @@
+###########################################################################
+#
+# 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.
+#
+###########################################################################
+if(iMSTK_USE_OpenHaptics)
+
+  project(Example-PBDCutting)
+
+  #-----------------------------------------------------------------------------
+  # Create executable
+  #-----------------------------------------------------------------------------
+  imstk_add_executable(${PROJECT_NAME} PBDCuttingExample.cpp)
+
+  #-----------------------------------------------------------------------------
+  # Add the target to Examples folder
+  #-----------------------------------------------------------------------------
+  SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples/Haptics)
+  
+  #-----------------------------------------------------------------------------
+  # Link libraries to executable
+  #-----------------------------------------------------------------------------
+  target_link_libraries(${PROJECT_NAME} SimulationManager apiUtilities Filtering)
+
+endif()
diff --git a/Examples/PBDCutting/PBDCuttingExample.cpp b/Examples/PBDCutting/PBDCuttingExample.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4048483ee8aa5acbc1313d7c82482d975061ea66
--- /dev/null
+++ b/Examples/PBDCutting/PBDCuttingExample.cpp
@@ -0,0 +1,243 @@
+/*=========================================================================
+
+   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 "imstkHapticDeviceClient.h"
+#include "imstkHapticDeviceManager.h"
+#include "imstkKeyboardDeviceClient.h"
+#include "imstkKeyboardSceneControl.h"
+#include "imstkLight.h"
+#include "imstkLogger.h"
+#include "imstkMouseSceneControl.h"
+#include "imstkNew.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkPbdObjectCuttingPair.h"
+#include "imstkRenderMaterial.h"
+#include "imstkScene.h"
+#include "imstkSceneManager.h"
+#include "imstkSceneObjectController.h"
+#include "imstkSimulationManager.h"
+#include "imstkSurfaceMesh.h"
+#include "imstkVisualModel.h"
+#include "imstkVTKViewer.h"
+
+using namespace imstk;
+
+// Parameters to play with
+const double width  = 50.0;
+const double height = 50.0;
+const int    nRows  = 12;
+const int    nCols  = 12;
+
+///
+/// \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("Cloth_SurfaceMesh");
+
+    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 [/][\] pattern
+            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.0e3);
+    pbdParams->enableConstraint(PbdConstraint::Type::Dihedral, 1.0e3);
+    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->setCollidingGeometry(clothMesh);
+    clothObj->setDynamicalModel(pbdModel);
+
+    return clothObj;
+}
+
+///
+/// \brief This example demonstrates the concept of PBD cutting
+/// for haptic interaction. NOTE: Requires GeoMagic Touch device
+///
+int
+main()
+{
+    // Setup logger (write to file and stdout)
+    Logger::startLogger();
+
+    // Scene
+    imstkNew<Scene> scene("PBDCutting");
+
+    // Create a cutting plane object in the scene
+    std::shared_ptr<SurfaceMesh> cutGeom(makeClothGeometry(40, 40, 2, 2));
+    cutGeom->setTranslation(Vec3d(-10, -20, 0));
+    cutGeom->updatePostTransformData();
+    imstkNew<CollidingObject> cutObj("CuttingObject");
+    cutObj->setVisualGeometry(cutGeom);
+    cutObj->setCollidingGeometry(cutGeom);
+    cutObj->getVisualModel(0)->getRenderMaterial()->setDisplayMode(RenderMaterial::DisplayMode::WireframeSurface);
+    scene->addSceneObject(cutObj);
+
+    // Create a pbd cloth object in the scene
+    std::shared_ptr<PbdObject> clothObj = makeClothObj("Cloth", width, height, nRows, nCols);
+    scene->addSceneObject(clothObj);
+
+    // Add interaction pair for pbd cutting
+    imstkNew<PbdObjectCuttingPair> cuttingPair(clothObj, cutObj);
+
+    // Device Server
+    imstkNew<HapticDeviceManager>       server;
+    std::shared_ptr<HapticDeviceClient> client = server->makeDeviceClient();
+
+    // Create the virtual coupling object controller
+    imstkNew<SceneObjectController> controller(cutObj, client);
+    scene->addController(controller);
+
+    // Camera
+    scene->getActiveCamera()->setPosition(Vec3d(1.0, 1.0, 1.0) * 100.0);
+    scene->getActiveCamera()->setFocalPoint(Vec3d(0, -50, 0));
+
+    // Light
+    imstkNew<DirectionalLight> light("light");
+    light->setFocalPoint(Vec3d(5.0, -8.0, -5.0));
+    light->setIntensity(1.0);
+    scene->addLight(light);
+
+    //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);
+        sceneManager->setExecutionType(Module::ExecutionType::ADAPTIVE);
+        sceneManager->pause(); // Start simulation paused
+
+        imstkNew<SimulationManager> driver;
+        driver->addModule(server);
+        driver->addModule(viewer);
+        driver->addModule(sceneManager);
+        driver->setDesiredDt(0.001);
+
+        // 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->setModuleDriver(driver);
+            viewer->addControl(keyControl);
+        }
+
+        // Queue keypress to be called after scene thread
+        queueConnect<KeyEvent>(viewer->getKeyboardDevice(), &KeyboardDeviceClient::keyPress, sceneManager, [&](KeyEvent* e)
+        {
+            // When i is pressed replace the PBD cloth with a cut one
+            if (e->m_key == 'i' && e->m_keyPressType == KEY_PRESS)
+            {
+                cuttingPair->apply();
+            }
+            });
+
+        driver->start();
+    }
+    return 0;
+}
diff --git a/Examples/PBDPicking/PBDPickingExample.cpp b/Examples/PBDPicking/PBDPickingExample.cpp
index d69e2f1f031ac086df791406fa892250b4cdd465..7fa7b148c83afa317e0fdea8f46bd192d6042482 100644
--- a/Examples/PBDPicking/PBDPickingExample.cpp
+++ b/Examples/PBDPicking/PBDPickingExample.cpp
@@ -158,7 +158,6 @@ main()
 
     // Scene
     imstkNew<Scene> scene("PBDPicking");
-    scene->getConfig()->writeTaskGraph = true;
 
     // Create the virtual coupling object controller
 
@@ -177,7 +176,7 @@ main()
     geomShaft->setOrientationAxis(Vec3d(0.0, 0.0, 1.0));
     geomShaft->setTranslation(Vec3d(0.0, 0.0, 10.0));
     imstkNew<CollidingObject> objShaft("ShaftObject");
-    objShaft->setVisualGeometry(geomShaft);
+    objShaft->setVisualGeometry(pivotSurfMesh);
     objShaft->setCollidingGeometry(geomShaft);
     scene->addSceneObject(objShaft);
 
@@ -187,7 +186,7 @@ main()
     geomUpperJaw->setRadius(2.0);
     geomUpperJaw->setOrientationAxis(Vec3d(0.0, 0.0, 1.0));
     imstkNew<CollidingObject> objUpperJaw("UpperJawObject");
-    objUpperJaw->setVisualGeometry(geomUpperJaw);
+    objUpperJaw->setVisualGeometry(upperSurfMesh);
     objUpperJaw->setCollidingGeometry(geomUpperJaw);
     scene->addSceneObject(objUpperJaw);
 
@@ -197,7 +196,7 @@ main()
     geomLowerJaw->setRadius(2.0);
     geomLowerJaw->setOrientationAxis(Vec3d(0.0, 0.0, 1.0));
     imstkNew<CollidingObject> objLowerJaw("LowerJawObject");
-    objLowerJaw->setVisualGeometry(geomLowerJaw);
+    objLowerJaw->setVisualGeometry(lowerSurfMesh);
     objLowerJaw->setCollidingGeometry(geomLowerJaw);
     scene->addSceneObject(objLowerJaw);
 
diff --git a/Examples/SDFHaptics/SDFHapticsExample.cpp b/Examples/SDFHaptics/SDFHapticsExample.cpp
index fece4a25e59f9be78d8b8af12ad4df2f5f35ffd4..60badea5d7b914d41f08bb8993f6faabd6aecc3f 100644
--- a/Examples/SDFHaptics/SDFHapticsExample.cpp
+++ b/Examples/SDFHaptics/SDFHapticsExample.cpp
@@ -99,6 +99,7 @@ main()
         // Setup a scene manager to advance the scene in its own thread
         imstkNew<SceneManager> sceneManager("Scene Manager");
         sceneManager->setActiveScene(scene);
+        sceneManager->setExecutionType(Module::ExecutionType::ADAPTIVE);
 
         imstkNew<SimulationManager> driver;
         driver->addModule(viewer);
diff --git a/Source/Controllers/imstkLaparoscopicToolController.cpp b/Source/Controllers/imstkLaparoscopicToolController.cpp
index 3d281ee99cae48494383e4d38e10af1d5ffa00c5..92c711c1dc73f87869d8e2af8927904b551b5b82 100644
--- a/Source/Controllers/imstkLaparoscopicToolController.cpp
+++ b/Source/Controllers/imstkLaparoscopicToolController.cpp
@@ -91,13 +91,19 @@ LaparoscopicToolController::update(const double dt)
 
     // TRS decompose and set upper/lower jaw geometries
     {
-        m_upperJaw->getVisualGeometry()->setTransform(m_controllerWorldTransform * m_upperJawLocalTransform * m_upperJawVisualTransform);
-        m_upperJaw->getCollidingGeometry()->setTransform(m_controllerWorldTransform * m_upperJawLocalTransform * m_upperJawCollidingTransform);
+        const Mat4d upperWorldTransform = m_controllerWorldTransform * m_upperJawLocalTransform;
+        m_upperJaw->getVisualGeometry()->setTransform(upperWorldTransform * m_upperJawVisualTransform);
+        m_upperJaw->getCollidingGeometry()->setTransform(upperWorldTransform * m_upperJawCollidingTransform);
     }
     {
-        m_lowerJaw->getVisualGeometry()->setTransform(m_controllerWorldTransform * m_lowerJawLocalTransform * m_upperJawVisualTransform);
-        m_lowerJaw->getCollidingGeometry()->setTransform(m_controllerWorldTransform * m_lowerJawLocalTransform * m_lowerJawCollidingTransform);
+        const Mat4d lowerWorldTransform = m_controllerWorldTransform * m_lowerJawLocalTransform;
+        m_lowerJaw->getVisualGeometry()->setTransform(lowerWorldTransform * m_upperJawVisualTransform);
+        m_lowerJaw->getCollidingGeometry()->setTransform(lowerWorldTransform * m_lowerJawCollidingTransform);
     }
+
+    m_shaft->getVisualGeometry()->postModified();
+    m_lowerJaw->getVisualGeometry()->postModified();
+    m_upperJaw->getVisualGeometry()->postModified();
 }
 
 void
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index b648246710bb2131c7fad93dba89a14f2be8305a..fd14a7c8eb8dd0ea0b1307ff79e05118d3484b93 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -36,6 +36,9 @@
 #include "imstkTaskGraph.h"
 #include "imstkTetrahedralMesh.h"
 
+#include <map>
+#include <set>
+
 namespace imstk
 {
 PbdModel::PbdModel() : DynamicalModel(DynamicalModelType::PositionBasedDynamics),
@@ -520,18 +523,24 @@ PbdModel::initializeDihedralConstraints(const double stiffness)
                 rs.resize(static_cast<size_t>(it - rs.begin()));
                 if (rs.size() > 1)
                 {
-                    size_t      idx = (rs[0] == k) ? 1 : 0;
-                    const auto& tri = elements[rs[idx]];
+                    size_t      idx  = (rs[0] == k) ? 1 : 0;
+                    const auto& tri0 = elements[k];
+                    const auto& tri1 = elements[rs[idx]];
+                    size_t      idx0 = 0;
+                    size_t      idx1 = 0;
                     for (size_t i = 0; i < 3; ++i)
                     {
-                        if (tri[i] != tri[0] && tri[i] != tri[1])
+                        if (tri0[i] != i1 && tri0[i] != i2)
                         {
-                            idx = i;
-                            break;
+                            idx0 = tri0[i];
+                        }
+                        if (tri1[i] != i1 && tri1[i] != i2)
+                        {
+                            idx1 = tri1[i];
                         }
                     }
                     auto c = std::make_shared<PbdDihedralConstraint>();
-                    c->initConstraint(*m_initialState->getPositions(), tri[2], tri[idx], tri[0], tri[1], stiffness);
+                    c->initConstraint(*m_initialState->getPositions(), idx0, idx1, i1, i2, stiffness);
                     m_constraints->push_back(std::move(c));
                 }
             }
@@ -570,6 +579,199 @@ PbdModel::initializeConstantDensityConstraint(const double stiffness)
     return true;
 }
 
+void
+PbdModel::removeConstraints(std::shared_ptr<std::unordered_set<size_t>> vertices)
+{
+    // constraint removal predicate
+    auto removeConstraint =
+        [&](std::shared_ptr<PbdConstraint> constraint)
+        {
+            for (auto i : constraint->getVertexIds())
+            {
+                if (vertices->find(i) != vertices->end())
+                {
+                    return true;
+                }
+            }
+            return false;
+        };
+
+    // remove constraints from constraint-partition maps and constraint vectors.
+    m_constraints->erase(std::remove_if(m_constraints->begin(), m_constraints->end(), removeConstraint),
+                         m_constraints->end());
+    for (auto& pc : *m_partitionedConstraints)
+    {
+        pc.erase(std::remove_if(pc.begin(), pc.end(), removeConstraint), pc.end());
+    }
+}
+
+void
+PbdModel::addConstraints(std::shared_ptr<std::unordered_set<size_t>> vertices)
+{
+    // check if constraint type matches the mesh type
+    CHECK(m_mesh->getTypeName() == "SurfaceMesh")
+        << "Add element constraints does not support current mesh type.";
+
+    const auto&                      triMesh  = std::static_pointer_cast<SurfaceMesh>(m_mesh);
+    const auto                       nV       = triMesh->getNumVertices();
+    const auto                       elements = triMesh->getTriangleIndices();
+    std::vector<std::vector<size_t>> onering(nV);
+
+    // build onering
+    for (auto& vertOnering : onering)
+    {
+        vertOnering.reserve(10);
+    }
+    for (int k = 0; k < elements->size(); ++k)
+    {
+        auto& tri = (*elements)[k];
+        onering[tri[0]].push_back(k);
+        onering[tri[1]].push_back(k);
+        onering[tri[2]].push_back(k);
+    }
+    for (auto& vertOnering : onering)
+    {
+        std::sort(vertOnering.begin(), vertOnering.end());
+    }
+
+    // functions for adding constraints
+    auto addDistanceConstraint =
+        [&](size_t i1, size_t i2, double stiffness)
+        {
+            auto c = std::make_shared<PbdDistanceConstraint>();
+            c->initConstraint(*m_initialState->getPositions(), i1, i2, stiffness);
+            m_constraints->push_back(std::move(c));
+        };
+    auto addAreaConstraint =
+        [&](size_t k, double stiffness)
+        {
+            auto& tri = (*elements)[k];
+            auto  c   = std::make_shared<PbdAreaConstraint>();
+            c->initConstraint(*m_initialState->getPositions(), tri[0], tri[1], tri[2], stiffness);
+            m_constraints->push_back(std::move(c));
+        };
+    auto addDihedralConstraint =
+        [&](size_t t0, size_t t1, size_t i1, size_t i2, double stiffness)
+        {
+            const auto& tri0 = (*elements)[t0];
+            const auto& tri1 = (*elements)[t1];
+            size_t      idx0 = 0;
+            size_t      idx1 = 0;
+            for (size_t i = 0; i < 3; ++i)
+            {
+                if (tri0[i] != i1 && tri0[i] != i2)
+                {
+                    idx0 = tri0[i];
+                }
+                if (tri1[i] != i1 && tri1[i] != i2)
+                {
+                    idx1 = tri1[i];
+                }
+            }
+            auto c = std::make_shared<PbdDihedralConstraint>();
+            c->initConstraint(*m_initialState->getPositions(), idx0, idx1, i1, i2, stiffness);
+            m_constraints->push_back(std::move(c));
+        };
+
+    // count constraints to be added for pre-allocation
+    std::set<std::pair<size_t, size_t>>                            distanceSet;
+    std::unordered_set<size_t>                                     areaSet;
+    std::map<std::pair<size_t, size_t>, std::pair<size_t, size_t>> dihedralSet;
+    for (auto& constraint : m_parameters->m_regularConstraints)
+    {
+        switch (constraint.first)
+        {
+        case PbdConstraint::Type::Distance:
+            for (const auto& vertIdx : *vertices)
+            {
+                for (const auto& triIdx : onering[vertIdx])
+                {
+                    const auto& tri = (*elements)[triIdx];
+                    size_t      i1  = 0;
+                    size_t      i2  = 0;
+                    for (size_t i = 0; i < 3; i++)
+                    {
+                        if (tri[i] == vertIdx)
+                        {
+                            i1 = tri[(i + 1) % 3];
+                            i2 = tri[(i + 2) % 3];
+                            break;
+                        }
+                    }
+                    auto pair1 = std::make_pair(std::min(vertIdx, i1), std::max(vertIdx, i1));
+                    auto pair2 = std::make_pair(std::min(vertIdx, i2), std::max(vertIdx, i2));
+                    distanceSet.insert(pair1);
+                    distanceSet.insert(pair2);
+                }
+            }
+            break;
+        case PbdConstraint::Type::Area:
+            for (const auto& vertIdx : *vertices)
+            {
+                for (const auto& triIdx : onering[vertIdx])
+                {
+                    areaSet.insert(triIdx);
+                }
+            }
+            break;
+        case PbdConstraint::Type::Dihedral:
+            for (const auto& vertIdx : *vertices)
+            {
+                for (const auto& triIdx : onering[vertIdx])
+                {
+                    const auto& tri = (*elements)[triIdx];
+                    for (size_t i = 0; i < 3; i++)
+                    {
+                        size_t j  = (i + 1) % 3;
+                        size_t i0 = tri[i];
+                        size_t i1 = tri[j];
+                        if (i0 > i1)
+                        {
+                            std::swap(i0, i1);
+                        }
+                        auto&               r0 = onering[i0];
+                        auto&               r1 = onering[i1];
+                        std::vector<size_t> rs(2);
+                        auto                it = std::set_intersection(r0.begin(), r0.end(), r1.begin(), r1.end(), rs.begin());
+                        rs.resize(static_cast<size_t>(it - rs.begin()));
+                        if (rs.size() > 1)
+                        {
+                            dihedralSet[std::make_pair(i0, i1)] = std::make_pair(rs[0], rs[1]);
+                        }
+                    }
+                }
+            }
+            break;
+        default:
+            break;
+        }
+    }
+
+    // add constraints
+    m_constraints->reserve(m_constraints->size() + distanceSet.size() + areaSet.size() + dihedralSet.size());
+    for (auto& constraint : m_parameters->m_regularConstraints)
+    {
+        switch (constraint.first)
+        {
+        case PbdConstraint::Type::Distance:
+            for (auto& c : distanceSet)
+            {
+                addDistanceConstraint(c.first, c.second, constraint.second);
+            }
+        case PbdConstraint::Type::Area:
+            for (auto& c : areaSet)
+            {
+                addAreaConstraint(c, constraint.second);
+            }
+        case PbdConstraint::Type::Dihedral:
+            for (auto& c : dihedralSet)
+            {
+                addDihedralConstraint(c.second.first, c.second.second, c.first.first, c.first.second, constraint.second);
+            }
+        }
+    }
+}
+
 void
 PbdModel::partitionConstraints(const bool print)
 {
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
index 71442b3232d8bf1d5f088b2efc082404118dadf9..e01e211c05c18ec0dce9156c3eb4761544c7489c 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -27,6 +27,7 @@
 #include "imstkPbdState.h"
 
 #include <unordered_map>
+#include <unordered_set>
 
 namespace imstk
 {
@@ -163,6 +164,17 @@ public:
     ///
     bool hasConstraints() const { return !m_constraints->empty() || !m_partitionedConstraints->empty(); }
 
+    ///
+    /// \brief Remove constraints related to a set of vertices.
+    ///
+    void removeConstraints(std::shared_ptr<std::unordered_set<size_t>> vertices);
+
+    ///
+    /// \brief Add constraints related to a set of vertices.
+    /// \brief Does not check for duplicating pre-existed constraints.
+    ///
+    void addConstraints(std::shared_ptr<std::unordered_set<size_t>> vertices);
+
     ///
     /// \brief Set the time step size
     ///
diff --git a/Source/Filtering/imstkSurfaceMeshCut.cpp b/Source/Filtering/imstkSurfaceMeshCut.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11c0f002267d19d7207b04b9447716d21c7b6e2c
--- /dev/null
+++ b/Source/Filtering/imstkSurfaceMeshCut.cpp
@@ -0,0 +1,668 @@
+/*=========================================================================
+
+   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 "imstkSurfaceMeshCut.h"
+
+#include "imstkGeometryUtilities.h"
+#include "imstkLogger.h"
+#include "imstkPlane.h"
+#include "imstkAnalyticalGeometry.h"
+
+#include <set>
+
+namespace imstk
+{
+SurfaceMeshCut::SurfaceMeshCut()
+{
+    setNumberOfInputPorts(1);
+    setNumberOfOutputPorts(1);
+    setOutput(std::make_shared<SurfaceMesh>());
+
+    m_CutGeometry = std::make_shared<Plane>();
+    m_RemoveConstraintVertices = std::make_shared<std::unordered_set<size_t>>();
+    m_AddConstraintVertices    = std::make_shared<std::unordered_set<size_t>>();
+}
+
+std::shared_ptr<SurfaceMesh>
+SurfaceMeshCut::getOutputMesh()
+{
+    return std::dynamic_pointer_cast<SurfaceMesh>(getOutput(0));
+}
+
+void
+SurfaceMeshCut::setInputMesh(std::shared_ptr<SurfaceMesh> inputMesh)
+{
+    setInput(inputMesh, 0);
+}
+
+void
+SurfaceMeshCut::requestUpdate()
+{
+    // input and output SurfaceMesh
+    std::shared_ptr<SurfaceMesh> inputSurf = std::dynamic_pointer_cast<SurfaceMesh>(getInput(0));
+    if (inputSurf == nullptr)
+    {
+        LOG(WARNING) << "Missing required SurfaceMesh input";
+        return;
+    }
+    std::shared_ptr<SurfaceMesh> outputSurf = std::dynamic_pointer_cast<SurfaceMesh>(getOutput(0));
+    outputSurf->deepCopy(inputSurf);
+
+    // vertices on the cutting path and whether they will be split
+    std::map<int, bool> cutVerts;
+
+    // generate cut data
+    if (std::dynamic_pointer_cast<AnalyticalGeometry>(m_CutGeometry) != nullptr)
+    {
+        generateAnalyticalCutData(std::static_pointer_cast<AnalyticalGeometry>(m_CutGeometry), outputSurf);
+    }
+    else if (std::dynamic_pointer_cast<SurfaceMesh>(m_CutGeometry) != nullptr)
+    {
+        generateSurfaceMeshCutData(std::static_pointer_cast<SurfaceMesh>(m_CutGeometry), outputSurf);
+    }
+    else
+    {
+        setOutput(outputSurf);
+        return;
+    }
+
+    if (m_CutData->size() == 0)
+    {
+        setOutput(outputSurf);
+        return;
+    }
+
+    // refinement
+    refinement(outputSurf, cutVerts);
+
+    // split cutting vertices
+    splitVerts(outputSurf, cutVerts, m_CutGeometry);
+
+    // vtkPolyData to output SurfaceMesh
+    setOutput(outputSurf);
+}
+
+void
+SurfaceMeshCut::refinement(std::shared_ptr<SurfaceMesh> outputSurf, std::map<int, bool>& cutVerts)
+{
+    // map between one exsiting edge to the new vert generated from the cutting
+    std::map<std::pair<int, int>, int> edgeVertMap;
+
+    auto triangles = outputSurf->getTriangleIndices();
+    auto vertices  = outputSurf->getVertexPositions();
+    auto initVerts = outputSurf->getInitialVertexPositions();
+    triangles->reserve(triangles->size() * 2);
+    vertices->reserve(vertices->size() * 2);
+    initVerts->reserve(initVerts->size() * 2);
+
+    for (const auto& curCutData : *m_CutData)
+    {
+        auto  curCutType = curCutData.cutType;
+        int   triId      = curCutData.triId;
+        int   ptId0      = curCutData.ptIds[0];
+        int   ptId1      = curCutData.ptIds[1];
+        Vec3d cood0      = curCutData.cutCoords[0];
+        Vec3d cood1      = curCutData.cutCoords[1];
+        Vec3d initCood0  = curCutData.initCoords[0];
+        Vec3d initCood1  = curCutData.initCoords[1];
+
+        if (curCutType == CutType::EDGE || curCutType == CutType::EDGE_VERT)
+        {
+            int  newPtId = -1;
+            auto edge0   = std::make_pair(ptId0, ptId1);
+            auto edge1   = std::make_pair(ptId1, ptId0);
+
+            // add new point
+            if (edgeVertMap.find(edge1) != edgeVertMap.end())
+            {
+                newPtId = edgeVertMap[edge1];
+            }
+            else
+            {
+                newPtId = vertices->size();
+                vertices->push_back(cood0);
+                initVerts->push_back(initCood0);
+                edgeVertMap[edge0] = newPtId;
+            }
+
+            // update triangle indices
+            Vec3i ptIds = (*triangles)[triId];
+            int   ptId2 = -1;
+            if (ptIds[0] != ptId0 && ptIds[0] != ptId1)
+            {
+                ptId2 = ptIds[0];
+            }
+            else if (ptIds[1] != ptId0 && ptIds[1] != ptId1)
+            {
+                ptId2 = ptIds[1];
+            }
+            else
+            {
+                ptId2 = ptIds[2];
+            }
+            (*triangles)[triId] = Vec3i(ptId2, ptId0, newPtId);
+            triangles->push_back(Vec3i(ptId2, newPtId, ptId1));
+
+            // add vertices to cutting path
+            if (curCutType == CutType::EDGE_VERT)
+            {
+                cutVerts[ptId2]   = (cutVerts.find(ptId2) != cutVerts.end()) ? true : false;
+                cutVerts[newPtId] = (cutVerts.find(newPtId) != cutVerts.end()) ? true : false;
+            }
+
+            m_RemoveConstraintVertices->insert(ptId0);
+            m_RemoveConstraintVertices->insert(ptId1);
+            m_RemoveConstraintVertices->insert(ptId2);
+            m_AddConstraintVertices->insert(ptId0);
+            m_AddConstraintVertices->insert(ptId1);
+            m_AddConstraintVertices->insert(ptId2);
+            m_AddConstraintVertices->insert(newPtId);
+        }
+        else if (curCutType == CutType::EDGE_EDGE)
+        {
+            int newPtId0 = -1;
+            int newPtId1 = -1;
+
+            Vec3i ptIds = (*triangles)[triId];
+            int   ptId2 = -1;
+            if (ptIds[0] != ptId0 && ptIds[0] != ptId1)
+            {
+                ptId2 = ptIds[0];
+            }
+            else if (ptIds[1] != ptId0 && ptIds[1] != ptId1)
+            {
+                ptId2 = ptIds[1];
+            }
+            else
+            {
+                ptId2 = ptIds[2];
+            }
+
+            auto edge00 = std::make_pair(ptId2, ptId0);
+            auto edge01 = std::make_pair(ptId0, ptId2);
+            auto edge10 = std::make_pair(ptId1, ptId2);
+            auto edge11 = std::make_pair(ptId2, ptId1);
+
+            // add new points
+            if (edgeVertMap.find(edge01) != edgeVertMap.end())
+            {
+                newPtId0 = edgeVertMap[edge01];
+            }
+            else
+            {
+                newPtId0 = vertices->size();
+                vertices->push_back(cood0);
+                initVerts->push_back(initCood0);
+                edgeVertMap[edge00] = newPtId0;
+            }
+            if (edgeVertMap.find(edge11) != edgeVertMap.end())
+            {
+                newPtId1 = edgeVertMap[edge11];
+            }
+            else
+            {
+                newPtId1 = vertices->size();
+                vertices->push_back(cood1);
+                initVerts->push_back(initCood1);
+                edgeVertMap[edge10] = newPtId1;
+            }
+
+            // update triangle indices
+            (*triangles)[triId] = Vec3i(ptId2, newPtId0, newPtId1);
+            triangles->push_back(Vec3i(newPtId0, ptId0, ptId1));
+            triangles->push_back(Vec3i(newPtId0, ptId1, newPtId1));
+
+            // add vertices to cutting path
+            cutVerts[newPtId0] = (cutVerts.find(newPtId0) != cutVerts.end()) ? true : false;
+            cutVerts[newPtId1] = (cutVerts.find(newPtId1) != cutVerts.end()) ? true : false;
+
+            m_RemoveConstraintVertices->insert(ptId0);
+            m_RemoveConstraintVertices->insert(ptId1);
+            m_RemoveConstraintVertices->insert(ptId2);
+            m_AddConstraintVertices->insert(ptId0);
+            m_AddConstraintVertices->insert(ptId1);
+            m_AddConstraintVertices->insert(ptId2);
+            m_AddConstraintVertices->insert(newPtId0);
+            m_AddConstraintVertices->insert(newPtId1);
+        }
+        else if (curCutType == CutType::VERT_VERT)
+        {
+            // add vertices to cutting path
+            cutVerts[ptId0] = (cutVerts.find(ptId0) != cutVerts.end()) ? true : false;
+            cutVerts[ptId1] = (cutVerts.find(ptId1) != cutVerts.end()) ? true : false;
+
+            m_RemoveConstraintVertices->insert(ptId0);
+            m_RemoveConstraintVertices->insert(ptId1);
+            m_AddConstraintVertices->insert(ptId0);
+            m_AddConstraintVertices->insert(ptId1);
+        }
+        else
+        {
+            //do nothing
+        }
+    }
+
+    outputSurf->modified();
+}
+
+void
+SurfaceMeshCut::splitVerts(std::shared_ptr<SurfaceMesh> outputSurf, std::map<int, bool>& cutVerts, std::shared_ptr<Geometry> geometry)
+{
+    auto triangles = outputSurf->getTriangleIndices();
+    auto vertices  = outputSurf->getVertexPositions();
+    auto initVerts = outputSurf->getInitialVertexPositions();
+
+    std::shared_ptr<AnalyticalGeometry> cutGeometry;
+    if (std::dynamic_pointer_cast<AnalyticalGeometry>(geometry) != nullptr)
+    {
+        cutGeometry = std::static_pointer_cast<AnalyticalGeometry>(geometry);
+    }
+    else if (std::dynamic_pointer_cast<SurfaceMesh>(geometry) != nullptr)
+    {
+        // assuming triangles in cutSurf are co-planar
+        auto cutSurf      = std::static_pointer_cast<SurfaceMesh>(geometry);
+        auto cutTriangles = cutSurf->getTriangleIndices();
+        auto cutVertices  = cutSurf->getVertexPositions();
+
+        // compute cutting plane (assuming all triangles in cutSurf are co-planar)
+        Vec3d p0 = (*cutVertices)[(*cutTriangles)[0][0]];
+        Vec3d p1 = (*cutVertices)[(*cutTriangles)[0][1]];
+        Vec3d p2 = (*cutVertices)[(*cutTriangles)[0][2]];
+        Vec3d cutNormal = ((p1 - p0).cross(p2 - p0)).normalized();
+        auto  cutPlane  = std::make_shared<Plane>(p0, cutNormal, "cutPlane");
+        cutGeometry = std::static_pointer_cast<AnalyticalGeometry>(cutPlane);
+    }
+
+    // build vertexToTriangleListMap
+    std::vector<std::set<int>> vertexNeighborTriangles;
+    vertexNeighborTriangles.clear();
+    vertexNeighborTriangles.resize(vertices->size());
+
+    int triangleId = 0;
+    for (const auto& tri : *triangles)
+    {
+        vertexNeighborTriangles.at(tri[0]).insert(triangleId);
+        vertexNeighborTriangles.at(tri[1]).insert(triangleId);
+        vertexNeighborTriangles.at(tri[2]).insert(triangleId);
+        triangleId++;
+    }
+
+    // split cutting vertices
+    for (const auto& cutVert : cutVerts)
+    {
+        if (cutVert.second == false && !vertexOnBoundary(triangles, vertexNeighborTriangles[cutVert.first]))
+        {
+            // do not split vertex since it's the cutting end in surface
+            (*m_CutVertMap)[cutVert.first] = cutVert.first;
+        }
+        else
+        {
+            // split vertex
+            int newPtId = vertices->size();
+            vertices->push_back((*vertices)[cutVert.first]);
+            initVerts->push_back((*initVerts)[cutVert.first]);
+            (*m_CutVertMap)[cutVert.first] = newPtId;
+            m_AddConstraintVertices->insert(newPtId);
+
+            for (const auto& t : vertexNeighborTriangles[cutVert.first])
+            {
+                //if triangle on the negative side
+                Vec3d pt0 = (*vertices)[(*triangles)[t][0]];
+                Vec3d pt1 = (*vertices)[(*triangles)[t][1]];
+                Vec3d pt2 = (*vertices)[(*triangles)[t][2]];
+
+                if (pointOnGeometrySide(pt0, cutGeometry) < 0
+                    || pointOnGeometrySide(pt1, cutGeometry) < 0
+                    || pointOnGeometrySide(pt2, cutGeometry) < 0)
+                {
+                    for (int i = 0; i < 3; i++)
+                    {
+                        if ((*triangles)[t][i] == cutVert.first)
+                        {
+                            (*triangles)[t][i] = newPtId;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    outputSurf->modified();
+}
+
+int
+SurfaceMeshCut::pointOnGeometrySide(Vec3d pt, std::shared_ptr<Geometry> geometry)
+{
+    if (std::dynamic_pointer_cast<AnalyticalGeometry>(geometry) != nullptr)
+    {
+        return pointOnAnalyticalSide(pt, std::static_pointer_cast<AnalyticalGeometry>(geometry));
+    }
+    else if (std::dynamic_pointer_cast<SurfaceMesh>(geometry) != nullptr)
+    {
+        // save for curve surface cutting
+    }
+    return 0;
+}
+
+int
+SurfaceMeshCut::pointOnAnalyticalSide(Vec3d pt, std::shared_ptr<AnalyticalGeometry> geometry)
+{
+    double normalProjection = geometry->getFunctionValue(pt);
+    if (normalProjection > m_Epsilon)
+    {
+        return 1;
+    }
+    else if (normalProjection < -m_Epsilon)
+    {
+        return -1;
+    }
+    else
+    {
+        return 0;
+    }
+    return 0;
+}
+
+bool
+SurfaceMeshCut::vertexOnBoundary(std::shared_ptr<VecDataArray<int, 3>> triangleIndices, std::set<int>& triSet)
+{
+    std::set<int> nonRepeatNeighborVerts;
+    for (const auto& tri : triSet)
+    {
+        for (int i = 0; i < 3; i++)
+        {
+            int ptId = (*triangleIndices)[tri][i];
+            if (nonRepeatNeighborVerts.find(ptId) != nonRepeatNeighborVerts.end())
+            {
+                nonRepeatNeighborVerts.erase(ptId);
+            }
+            else
+            {
+                nonRepeatNeighborVerts.insert(ptId);
+            }
+        }
+    }
+    return (nonRepeatNeighborVerts.size() >= 2 ? true : false);
+}
+
+void
+SurfaceMeshCut::generateAnalyticalCutData(std::shared_ptr<AnalyticalGeometry> geometry, std::shared_ptr<SurfaceMesh> outputSurf)
+{
+    auto triangles = outputSurf->getTriangleIndices();
+    auto vertices  = outputSurf->getVertexPositions();
+    auto initVerts = outputSurf->getInitialVertexPositions();
+
+    m_CutData->clear();
+    std::set<std::pair<int, int>> repeatEdges;  // make sure not boundary edge in vert-vert case
+
+    for (int i = 0; i < triangles->size(); i++)
+    {
+        auto&   tri = (*triangles)[i];
+        CutData newCutData;
+
+        auto ptSide = Vec3i(pointOnAnalyticalSide((*vertices)[tri[0]], geometry),
+                            pointOnAnalyticalSide((*vertices)[tri[1]], geometry),
+                            pointOnAnalyticalSide((*vertices)[tri[2]], geometry));
+
+        switch (ptSide.squaredNorm())
+        {
+        case 1:
+            for (int j = 0; j < 3; j++)
+            {
+                if (ptSide[j] != 0)
+                {
+                    int idx0  = (j + 1) % 3;
+                    int idx1  = (j + 2) % 3;
+                    int ptId0 = tri[idx0];
+                    int ptId1 = tri[idx1];
+
+                    std::pair<int, int> cutEdge = std::make_pair(ptId1, ptId0);
+                    // if find cut triangle from the other side of the edge, add newCutData
+                    if (repeatEdges.find(cutEdge) != repeatEdges.end())
+                    {
+                        newCutData.cutType       = CutType::VERT_VERT;
+                        newCutData.triId         = i;
+                        newCutData.ptIds[0]      = ptId0;
+                        newCutData.ptIds[1]      = ptId1;
+                        newCutData.cutCoords[0]  = (*vertices)[ptId0];
+                        newCutData.cutCoords[1]  = (*vertices)[ptId1];
+                        newCutData.initCoords[0] = (*initVerts)[ptId0];
+                        newCutData.initCoords[1] = (*initVerts)[ptId1];
+                        m_CutData->push_back(newCutData);
+                    }
+                    else
+                    {
+                        repeatEdges.insert(std::make_pair(ptId0, ptId1));
+                    }
+                }
+            }
+            break;
+        case 2:
+            if (ptSide.sum() == 0)
+            {
+                for (int j = 0; j < 3; j++)
+                {
+                    if (ptSide[j] == 0)
+                    {
+                        int    idx0     = (j + 1) % 3;
+                        int    idx1     = (j + 2) % 3;
+                        int    ptId0    = tri[idx0];
+                        int    ptId1    = tri[idx1];
+                        Vec3d  pos0     = (*vertices)[ptId0];
+                        Vec3d  pos1     = (*vertices)[ptId1];
+                        Vec3d  initPos0 = (*initVerts)[ptId0];
+                        Vec3d  initPos1 = (*initVerts)[ptId1];
+                        double func0    = geometry->getFunctionValue(pos0);
+                        double func1    = geometry->getFunctionValue(pos1);
+                        double frac     = -func0 / (func1 - func0);
+
+                        newCutData.cutType       = CutType::EDGE_VERT;
+                        newCutData.triId         = i;
+                        newCutData.ptIds[0]      = ptId0;
+                        newCutData.ptIds[1]      = ptId1;
+                        newCutData.cutCoords[0]  =  frac * (pos1 - pos0) + pos0;
+                        newCutData.cutCoords[1]  = (*vertices)[tri[j]];
+                        newCutData.initCoords[0] = frac * (initPos1 - initPos0) + initPos0;
+                        newCutData.initCoords[1] = (*initVerts)[tri[j]];
+                        m_CutData->push_back(newCutData);
+                    }
+                }
+            }
+            else
+            {
+                // newCutData.cutType = CutType::VERT;
+            }
+            break;
+        case 3:
+            if (ptSide.sum() == -1 || ptSide.sum() == 1)
+            {
+                for (int j = 0; j < 3; j++)
+                {
+                    if (ptSide[j] == -ptSide.sum())
+                    {
+                        int    idx0     = (j + 1) % 3;
+                        int    idx1     = (j + 2) % 3;
+                        int    ptId0    = tri[idx0];
+                        int    ptId1    = tri[idx1];
+                        int    ptId2    = tri[j];
+                        Vec3d  pos0     = (*vertices)[ptId0];
+                        Vec3d  pos1     = (*vertices)[ptId1];
+                        Vec3d  pos2     = (*vertices)[ptId2];
+                        Vec3d  initPos0 = (*initVerts)[ptId0];
+                        Vec3d  initPos1 = (*initVerts)[ptId1];
+                        Vec3d  initPos2 = (*initVerts)[ptId2];
+                        double func0    = geometry->getFunctionValue(pos0);
+                        double func1    = geometry->getFunctionValue(pos1);
+                        double func2    = geometry->getFunctionValue(pos2);
+                        double frac0    = -func0 / (func2 - func0);
+                        double frac1    = -func1 / (func2 - func1);
+
+                        newCutData.cutType       = CutType::EDGE_EDGE;
+                        newCutData.triId         = i;
+                        newCutData.ptIds[0]      = ptId0;
+                        newCutData.ptIds[1]      = ptId1;
+                        newCutData.cutCoords[0]  =  frac0 * (pos2 - pos0) + pos0;
+                        newCutData.cutCoords[1]  =  frac1 * (pos2 - pos1) + pos1;
+                        newCutData.initCoords[0] = frac0 * (initPos2 - initPos0) + initPos0;
+                        newCutData.initCoords[1] = frac1 * (initPos2 - initPos1) + initPos1;
+                        m_CutData->push_back(newCutData);
+                    }
+                }
+            }
+            else
+            {
+                // no intersection
+            }
+            break;
+        default:
+            break;
+        }
+    }
+}
+
+void
+SurfaceMeshCut::generateSurfaceMeshCutData(std::shared_ptr<SurfaceMesh> cutSurf, std::shared_ptr<SurfaceMesh> outputSurf)
+{
+    auto cutTriangles = cutSurf->getTriangleIndices();
+    auto cutVertices  = cutSurf->getVertexPositions();
+    auto triangles    = outputSurf->getTriangleIndices();
+    auto vertices     = outputSurf->getVertexPositions();
+
+    // compute cutting plane (assuming all triangles in cutSurf are co-planar)
+    Vec3d p0 = (*cutVertices)[(*cutTriangles)[0][0]];
+    Vec3d p1 = (*cutVertices)[(*cutTriangles)[0][1]];
+    Vec3d p2 = (*cutVertices)[(*cutTriangles)[0][2]];
+    Vec3d cutNormal = (p1 - p0).cross(p2 - p0);
+    auto  cutPlane  = std::make_shared<Plane>(p0, cutNormal, "cutPlane");
+
+    // Compute cut data using infinite cutPlane
+    generateAnalyticalCutData(std::static_pointer_cast<AnalyticalGeometry>(cutPlane), outputSurf);
+
+    // remove cutData that are out of the cutSurf
+    std::shared_ptr<std::vector<CutData>> newCutData = std::make_shared<std::vector<CutData>>();
+    for (const auto& curCutData : *m_CutData)
+    {
+        bool    coord0In      = pointProjectionInSurface(curCutData.cutCoords[0], cutSurf);
+        bool    coord1In      = pointProjectionInSurface(curCutData.cutCoords[1], cutSurf);
+        CutData newCurCutData = curCutData;
+
+        switch (curCutData.cutType)
+        {
+        case CutType::VERT_VERT:
+            if (coord0In && coord1In)
+            {
+                newCutData->push_back(newCurCutData);
+            }
+            break;
+        case CutType::EDGE_VERT:
+            if (coord0In)     // edge inside
+            {
+                if (coord1In) // vert inside
+                {
+                    newCutData->push_back(newCurCutData);
+                }
+                else
+                {
+                    newCurCutData.cutType = CutType::EDGE;
+                    newCutData->push_back(newCurCutData);
+                }
+            }
+            break;
+        case CutType::EDGE_EDGE:
+            if (coord0In)
+            {
+                if (coord1In)
+                {
+                    newCutData->push_back(newCurCutData);
+                }
+                else
+                {
+                    auto& tri = (*triangles)[newCurCutData.triId];
+                    for (int i = 0; i < 3; i++)
+                    {
+                        if (tri[i] == newCurCutData.ptIds[0])
+                        {
+                            int idx0 = (i + 2) % 3;
+                            newCurCutData.ptIds[0] = tri[idx0];
+                            newCurCutData.ptIds[1] = tri[i];
+                            break;
+                        }
+                    }
+                    newCurCutData.cutType = CutType::EDGE;
+                    newCutData->push_back(newCurCutData);
+                }
+            }
+            else if (coord1In) // && coord0Out
+            {
+                auto& tri = (*triangles)[newCurCutData.triId];
+                for (int i = 0; i < 3; i++)
+                {
+                    if (tri[i] == curCutData.ptIds[0])
+                    {
+                        int idx0 = (i + 1) % 3;
+                        int idx1 = (i + 2) % 3;
+                        newCurCutData.ptIds[0] = tri[idx0];
+                        newCurCutData.ptIds[1] = tri[idx1];
+                        break;
+                    }
+                }
+                newCurCutData.cutCoords[0]  = newCurCutData.cutCoords[1];
+                newCurCutData.initCoords[0] = newCurCutData.initCoords[1];
+                newCurCutData.cutType       = CutType::EDGE;
+                newCutData->push_back(newCurCutData);
+            }
+            break;
+        default:
+            break;
+        }
+    }
+
+    // update cutData
+    m_CutData = newCutData;
+}
+
+bool
+SurfaceMeshCut::pointProjectionInSurface(const Vec3d& pt, std::shared_ptr<SurfaceMesh> cutSurf)
+{
+    auto cutTriangles = cutSurf->getTriangleIndices();
+    auto cutVertices  = cutSurf->getVertexPositions();
+    bool inSurface    = false;
+
+    for (const auto& tri : *cutTriangles)
+    {
+        Vec3d p0     = (*cutVertices)[tri[0]];
+        Vec3d p1     = (*cutVertices)[tri[1]];
+        Vec3d p2     = (*cutVertices)[tri[2]];
+        Vec3d normal = (p1 - p0).cross(p2 - p0).normalized();
+
+        double leftP0P1 = normal.dot((p1 - p0).cross(pt - p0));
+        double leftP1P2 = normal.dot((p2 - p1).cross(pt - p1));
+        double leftP2P0 = normal.dot((p0 - p2).cross(pt - p2));
+
+        if (leftP0P1 >= 0 && leftP1P2 >= 0 && leftP2P0 >= 0)
+        {
+            inSurface = true;
+            break;
+        }
+    }
+    return inSurface;
+}
+}
diff --git a/Source/Filtering/imstkSurfaceMeshCut.h b/Source/Filtering/imstkSurfaceMeshCut.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a415280433809a51c1c1d5375ca7ee4e7061aef
--- /dev/null
+++ b/Source/Filtering/imstkSurfaceMeshCut.h
@@ -0,0 +1,151 @@
+/*=========================================================================
+
+   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.
+
+=========================================================================*/
+#pragma once
+
+#include "imstkGeometryAlgorithm.h"
+#include "imstkSurfaceMesh.h"
+#include <map>
+#include <unordered_set>
+
+namespace imstk
+{
+class Geometry;
+class AnalyticalGeometry;
+
+// vertex on the plane (0), positive side (+1), negative side (-1)
+// pt0 and pt1 follows the triangle's indexing order when tri is presented
+// c0 and c1 are cutting coordinates stored in cutData
+enum class CutType
+{
+    NONE = 0,
+    /* triangle is not cut through
+    *       pt0 (-+1)
+    *           /  \
+    *       c0 /    \
+    *         / tri  \
+    * pt1 (+-1)------(?)
+    */
+    EDGE,
+    /*
+    *      (-+1)
+    *       /  \
+    *      /    \
+    *     / tri  \
+    *  (-+1)------(0) pt0/c0
+    */
+    VERT,
+    /*
+    *        (+-1) pt1
+    *        /  \
+    *    c1 /    \
+    *      / tri  \
+    *  (-+1)--c0--(+-1) pt0
+    */
+    EDGE_EDGE,
+    /*
+    *        pt0 (+-1)
+    *            /  \
+    *        c0 /    \
+    *          / tri  \
+    *  pt1 (-+1)------(0) c1
+    */
+    EDGE_VERT,
+    /*
+    * pt0/c0 (0)------(+-1)
+    *        /  \      /
+    *       /    \    /
+    *      /      \  /
+    *   (-+1)------(0) pt1/c1
+    */
+    VERT_VERT
+};
+
+struct CutData
+{
+    public:
+        Vec3d cutCoords[2];
+        Vec3d initCoords[2];
+        int triId       = -1;
+        int ptIds[2]    = { -1, -1 };
+        CutType cutType = CutType::NONE;
+};
+
+///
+/// \class SurfaceMeshCut
+///
+/// \brief This filter cuts the triangles of a SurfaceMesh into smaller
+/// triangles using input cutData
+/// \todo: test
+///
+class SurfaceMeshCut : public GeometryAlgorithm
+{
+public:
+    SurfaceMeshCut();
+    virtual ~SurfaceMeshCut() override = default;
+
+public:
+    std::shared_ptr<SurfaceMesh> getOutputMesh();
+    void setInputMesh(std::shared_ptr<SurfaceMesh> inputSurf);
+    std::shared_ptr<std::map<int, int>> getCutVertMap() { return m_CutVertMap; }
+
+    imstkGetMacro(CutData, std::shared_ptr<std::vector<CutData>>);
+    imstkSetMacro(CutData, std::shared_ptr<std::vector<CutData>>);
+    imstkGetMacro(CutGeometry, std::shared_ptr<Geometry>);
+    imstkSetMacro(CutGeometry, std::shared_ptr<Geometry>);
+    imstkGetMacro(Epsilon, double);
+    imstkSetMacro(Epsilon, double);
+    imstkGetMacro(RemoveConstraintVertices, std::shared_ptr<std::unordered_set<size_t>>);
+    imstkGetMacro(AddConstraintVertices, std::shared_ptr<std::unordered_set<size_t>>);
+
+protected:
+    void requestUpdate() override;
+
+    void refinement(std::shared_ptr<SurfaceMesh> outputSurf,
+                    std::map<int, bool>& cutVerts);
+
+    void splitVerts(std::shared_ptr<SurfaceMesh> outputSurf,
+                    std::map<int, bool>& cutVerts,
+                    std::shared_ptr<Geometry> geometry);
+
+    int pointOnGeometrySide(Vec3d pt, std::shared_ptr<Geometry> geometry);
+
+    int pointOnAnalyticalSide(Vec3d pt, std::shared_ptr<AnalyticalGeometry> geometry);
+
+    bool vertexOnBoundary(std::shared_ptr<VecDataArray<int, 3>> triangleIndices,
+                          std::set<int>& triSet);
+
+    bool pointProjectionInSurface(const Vec3d& pt, std::shared_ptr<SurfaceMesh> cutSurf);
+
+    void generateAnalyticalCutData(std::shared_ptr<AnalyticalGeometry> geometry,
+                                   std::shared_ptr<SurfaceMesh>        outputSurf);
+
+    void generateSurfaceMeshCutData(std::shared_ptr<SurfaceMesh> cutSurf,
+                                    std::shared_ptr<SurfaceMesh> outputSurf);
+
+private:
+    std::shared_ptr<std::vector<CutData>> m_CutData    = std::make_shared<std::vector<CutData>>();
+    std::shared_ptr<std::map<int, int>>   m_CutVertMap = std::make_shared<std::map<int, int>>();
+    std::shared_ptr<Geometry> m_CutGeometry = nullptr;
+    std::shared_ptr<std::unordered_set<size_t>> m_RemoveConstraintVertices = nullptr;
+    std::shared_ptr<std::unordered_set<size_t>> m_AddConstraintVertices    = nullptr;
+    double m_Epsilon = 1;
+};
+}
\ No newline at end of file
diff --git a/Source/Geometry/Analytic/imstkPlane.h b/Source/Geometry/Analytic/imstkPlane.h
index 129318900eb48c182fd5513e4031fd822e492e47..d0b6966e18bd0dbfd81a551d3b251ff6f0543d86 100644
--- a/Source/Geometry/Analytic/imstkPlane.h
+++ b/Source/Geometry/Analytic/imstkPlane.h
@@ -85,7 +85,7 @@ public:
     ///
     /// \brief Returns signed distance to surface at pos
     ///
-    double getFunctionValue(const Vec3d& pos) const override { return m_orientationAxis.dot(pos - m_position); }
+    double getFunctionValue(const Vec3d& pos) const override { return m_orientationAxisPostTransform.dot(pos - m_positionPostTransform); }
 
     ///
     /// \brief Get the min, max of the AABB around the plane
diff --git a/Source/Geometry/Implicit/imstkSignedDistanceField.cpp b/Source/Geometry/Implicit/imstkSignedDistanceField.cpp
index e32bb168f8c6f094c11f1d3e5c56fb0dd2259adc..2d149531f98a7bc6b5718be0f3e989ea9b784123 100644
--- a/Source/Geometry/Implicit/imstkSignedDistanceField.cpp
+++ b/Source/Geometry/Implicit/imstkSignedDistanceField.cpp
@@ -112,7 +112,8 @@ SignedDistanceField::getFunctionValue(const Vec3d& pos) const
     }
     else
     {
-        return std::numeric_limits<double>::min();
+        // If outside of the bounds, return positive (assume not inside)
+        return IMSTK_DOUBLE_MAX;
     }
 }
 
diff --git a/Source/Geometry/Mesh/imstkPointSet.cpp b/Source/Geometry/Mesh/imstkPointSet.cpp
index 3db88629081984d785f23e130651778a2a42ae59..a0f3926c8fc44d7e8324965ad07119fcab27d335 100644
--- a/Source/Geometry/Mesh/imstkPointSet.cpp
+++ b/Source/Geometry/Mesh/imstkPointSet.cpp
@@ -107,8 +107,8 @@ PointSet::getInitialVertexPosition(const size_t vertNum)
 void
 PointSet::setVertexPositions(std::shared_ptr<VecDataArray<double, 3>> vertices)
 {
-    m_vertexPositions  = vertices;
-    m_transformApplied = false;
+    m_vertexPositions = vertices;
+    //m_transformApplied = false;
 
     this->updatePostTransformData();
 }
diff --git a/Source/Scene/CMakeLists.txt b/Source/Scene/CMakeLists.txt
index 114253e60cc34d7e2dac0f43841a29db2d0920f5..5d1a4640b5f69e59067630d2063dab2976020dd0 100644
--- a/Source/Scene/CMakeLists.txt
+++ b/Source/Scene/CMakeLists.txt
@@ -9,6 +9,7 @@ imstk_add_library( Scene
     CollisionHandling
     SceneEntities
     Controllers
+    Filtering
   )
 
 #-----------------------------------------------------------------------------
diff --git a/Source/Scene/imstkObjectInteractionFactory.cpp b/Source/Scene/imstkObjectInteractionFactory.cpp
index 80b9cd1ece2734e9f346739e61373aea98e776ea..584dfc5c204706cd8993407c6ef067afde351f7b 100644
--- a/Source/Scene/imstkObjectInteractionFactory.cpp
+++ b/Source/Scene/imstkObjectInteractionFactory.cpp
@@ -26,6 +26,7 @@ limitations under the License.
 #include "imstkFeDeformableObject.h"
 #include "imstkPbdObject.h"
 #include "imstkPbdObjectCollisionPair.h"
+#include "imstkPbdObjectCuttingPair.h"
 #include "imstkPbdObjectPickingPair.h"
 #include "imstkPBDPickingCH.h"
 #include "imstkPenaltyCH.h"
@@ -53,9 +54,13 @@ makeObjectInteractionPair(std::shared_ptr<CollidingObject> obj1, std::shared_ptr
     {
         results = std::make_shared<PbdObjectCollisionPair>(std::dynamic_pointer_cast<PbdObject>(obj1), std::dynamic_pointer_cast<PbdObject>(obj2), cdType);
     }
+    else if (intType == InteractionType::PbdObjToCollidingObjCutting && isType<PbdObject>(obj1))
+    {
+        results = std::make_shared<PbdObjectCuttingPair>(std::dynamic_pointer_cast<PbdObject>(obj1), obj2);
+    }
     else if (intType == InteractionType::PbdObjToCollidingObjPicking && isType<PbdObject>(obj1))
     {
-        results = std::make_shared<PbdObjectPickingPair>(std::dynamic_pointer_cast<PbdObject>(obj1), std::dynamic_pointer_cast<CollidingObject>(obj2), cdType);
+        results = std::make_shared<PbdObjectPickingPair>(std::dynamic_pointer_cast<PbdObject>(obj1), obj2, cdType);
     }
     //else if (intType == InteractionType::PbdObjToCollidingObjCollision && isType<PbdObject>(obj1))
     //{
diff --git a/Source/Scene/imstkObjectInteractionFactory.h b/Source/Scene/imstkObjectInteractionFactory.h
index e24cb78f6ddbf11db61ad63c3450e39b72306a4c..09a3a403d06416673ea5974d8570988a9b7c777c 100644
--- a/Source/Scene/imstkObjectInteractionFactory.h
+++ b/Source/Scene/imstkObjectInteractionFactory.h
@@ -33,6 +33,7 @@ enum class InteractionType
 {
     PbdObjToPbdObjCollision,
     PbdObjToCollidingObjPicking,
+    PbdObjToCollidingObjCutting,
 
     PbdObjToCollidingObjCollision,
     SphObjToCollidingObjCollision,
diff --git a/Source/Scene/imstkPbdObjectCuttingPair.cpp b/Source/Scene/imstkPbdObjectCuttingPair.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..12b04c31f21ccff1856ea857ba604a28a7917dca
--- /dev/null
+++ b/Source/Scene/imstkPbdObjectCuttingPair.cpp
@@ -0,0 +1,190 @@
+/*=========================================================================
+
+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 "imstkPbdObjectCuttingPair.h"
+
+#include "imstkAnalyticalGeometry.h"
+#include "imstkCollidingObject.h"
+#include "imstkLogger.h"
+#include "imstkNew.h"
+#include "imstkPbdModel.h"
+#include "imstkPbdObject.h"
+#include "imstkPbdSolver.h"
+#include "imstkSurfaceMesh.h"
+#include "imstkSurfaceMeshCut.h"
+#include "imstkVecDataArray.h"
+
+namespace imstk
+{
+PbdObjectCuttingPair::PbdObjectCuttingPair(std::shared_ptr<PbdObject> pbdObj, std::shared_ptr<CollidingObject> cutObj) : ObjectInteractionPair(pbdObj, cutObj)
+{
+    // check whether pbd object is a surfacemesh
+    if (std::dynamic_pointer_cast<SurfaceMesh>(pbdObj->getPhysicsGeometry()) == nullptr)
+    {
+        LOG(WARNING) << "PbdObj is not a SurfaceMesh, could not create cutting pair";
+        return;
+    }
+
+    // check whether cut object is valid
+    if (std::dynamic_pointer_cast<SurfaceMesh>(cutObj->getCollidingGeometry()) == nullptr
+        && std::dynamic_pointer_cast<AnalyticalGeometry>(cutObj->getCollidingGeometry()) == nullptr)
+    {
+        LOG(WARNING) << "CutObj is neither a SurfaceMesh nor an AnalyticalGeometry, could not create cutting pair";
+        return;
+    }
+}
+
+void
+PbdObjectCuttingPair::apply()
+{
+    auto pbdObj   = std::static_pointer_cast<PbdObject>(m_objects.first);
+    auto cutObj   = std::static_pointer_cast<CollidingObject>(m_objects.second);
+    auto pbdModel = pbdObj->getPbdModel();
+    auto pbdMesh  = std::static_pointer_cast<SurfaceMesh>(pbdModel->getModelGeometry());
+
+    m_addConstraintVertices->clear();
+    m_removeConstraintVertices->clear();
+
+    // Perform cutting
+    imstkNew<SurfaceMeshCut> surfCut;
+    surfCut->setInputMesh(pbdMesh);
+    surfCut->setCutGeometry(cutObj->getCollidingGeometry());
+    surfCut->update();
+    auto newPbdMesh = surfCut->getOutputMesh();
+
+    // Only remove and add constraints related to the topological changes
+    m_removeConstraintVertices = surfCut->getRemoveConstraintVertices();
+    m_addConstraintVertices    = surfCut->getAddConstraintVertices();
+
+    // update pbd mesh
+    pbdMesh->setInitialVertexPositions(std::make_shared<VecDataArray<double, 3>>(*newPbdMesh->getInitialVertexPositions()));
+    pbdMesh->setVertexPositions(std::make_shared<VecDataArray<double, 3>>(*newPbdMesh->getVertexPositions()));
+    pbdMesh->setTriangleIndices(std::make_shared<VecDataArray<int, 3>>(*newPbdMesh->getTriangleIndices()));
+    pbdMesh->modified();
+
+    // update pbd states, constraints and solver
+    pbdModel->initState();
+    pbdModel->removeConstraints(m_removeConstraintVertices);
+    pbdModel->addConstraints(m_addConstraintVertices);
+    pbdModel->getSolver()->setInvMasses(pbdModel->getInvMasses());
+    pbdModel->getSolver()->setPositions(pbdModel->getCurrentState()->getPositions());
+}
+
+void
+PbdObjectCuttingPair::addVertices(std::shared_ptr<SurfaceMesh> pbdMesh,
+                                  std::shared_ptr<VecDataArray<double, 3>> newVertices,
+                                  std::shared_ptr<VecDataArray<double, 3>> newInitialVertices)
+{
+    auto vertices = pbdMesh->getVertexPositions();
+    auto initialVertices = pbdMesh->getInitialVertexPositions();
+
+    auto nVertices    = vertices->size();
+    auto nNewVertices = newVertices->size();
+    if (nNewVertices != newInitialVertices->size())
+    {
+        LOG(WARNING) << "Number of new vertices does not match number of new initial vertices";
+        return;
+    }
+
+    vertices->reserve(nVertices + nNewVertices);
+    initialVertices->reserve(nVertices + nNewVertices);
+    for (int i = 0; i < nNewVertices; ++i)
+    {
+        vertices->push_back((*newVertices)[i]);
+        initialVertices->push_back((*newInitialVertices)[i]);
+    }
+}
+
+void
+PbdObjectCuttingPair::modifyVertices(std::shared_ptr<SurfaceMesh> pbdMesh,
+                                     std::shared_ptr<std::vector<size_t>> modifiedVertexIndices,
+                                     std::shared_ptr<VecDataArray<double, 3>> modifiedVertices,
+                                     std::shared_ptr<VecDataArray<double, 3>> modifiedInitialVertices)
+{
+    auto vertices = pbdMesh->getVertexPositions();
+    auto initialVertices = pbdMesh->getInitialVertexPositions();
+
+    auto nModifiedVertices = modifiedVertices->size();
+    if (nModifiedVertices != modifiedInitialVertices->size()
+        || nModifiedVertices != modifiedVertexIndices->size())
+    {
+        LOG(WARNING) << "Numbers of vertices do not match.";
+        return;
+    }
+
+    for (int i = 0; i < nModifiedVertices; ++i)
+    {
+        auto vertexIdx = modifiedVertexIndices->at(i);
+        (*vertices)[vertexIdx] = (*modifiedVertices)[i];
+        (*initialVertices)[vertexIdx] = (*modifiedInitialVertices)[i];
+        m_removeConstraintVertices->insert(vertexIdx);
+        m_addConstraintVertices->insert(vertexIdx);
+    }
+}
+
+void
+PbdObjectCuttingPair::addTriangles(std::shared_ptr<SurfaceMesh> pbdMesh,
+                                   std::shared_ptr<VecDataArray<int, 3>> newTriangles)
+{
+    auto triangles     = pbdMesh->getTriangleIndices();
+    auto nTriangles    = triangles->size();
+    auto nNewTriangles = newTriangles->size();
+
+    triangles->reserve(nTriangles + nNewTriangles);
+    for (int i = 0; i < nNewTriangles; ++i)
+    {
+        auto& tri = (*newTriangles)[i];
+        triangles->push_back(tri);
+        m_addConstraintVertices->insert(tri[0]);
+        m_addConstraintVertices->insert(tri[1]);
+        m_addConstraintVertices->insert(tri[2]);
+    }
+}
+
+void
+PbdObjectCuttingPair::modifyTriangles(std::shared_ptr<SurfaceMesh> pbdMesh,
+                                      std::shared_ptr<std::vector<size_t>> modifiedTriangleIndices,
+                                      std::shared_ptr<VecDataArray<int, 3>> modifiedTriangles)
+{
+    auto triangles = pbdMesh->getTriangleIndices();
+    auto nModifiedTriangles = modifiedTriangles->size();
+    if (nModifiedTriangles != modifiedTriangleIndices->size())
+    {
+        LOG(WARNING) << "Numbers of vertices do not match.";
+        return;
+    }
+
+    for (int i = 0; i < nModifiedTriangles; ++i)
+    {
+        auto  triId  = (*modifiedTriangleIndices)[i];
+        auto& oldTri = (*triangles)[triId];
+        m_removeConstraintVertices->insert(oldTri[0]);
+        m_removeConstraintVertices->insert(oldTri[1]);
+        m_removeConstraintVertices->insert(oldTri[2]);
+
+        auto& newTri = (*modifiedTriangles)[i];
+        (*triangles)[triId] = newTri;
+        m_addConstraintVertices->insert(newTri[0]);
+        m_addConstraintVertices->insert(newTri[1]);
+        m_addConstraintVertices->insert(newTri[2]);
+    }
+}
+}
\ No newline at end of file
diff --git a/Source/Scene/imstkPbdObjectCuttingPair.h b/Source/Scene/imstkPbdObjectCuttingPair.h
new file mode 100644
index 0000000000000000000000000000000000000000..f045183f2ce511a857116b345e06f688155d3669
--- /dev/null
+++ b/Source/Scene/imstkPbdObjectCuttingPair.h
@@ -0,0 +1,83 @@
+/*=========================================================================
+
+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.
+
+=========================================================================*/
+
+#pragma once
+
+#include "imstkObjectInteractionPair.h"
+#include "imstkVecDataArray.h"
+
+#include <unordered_set>
+#include <vector>
+
+namespace imstk
+{
+class CollidingObject;
+class PbdObject;
+class SurfaceMesh;
+
+///
+/// \class PbdObjectCuttingPair
+///
+/// \brief This class defines a cutting interaction between a PbdObject and a CollidingObject
+/// call apply to perform the cut given the current states of both objects
+///
+class PbdObjectCuttingPair : public ObjectInteractionPair
+{
+public:
+    PbdObjectCuttingPair(std::shared_ptr<PbdObject> pbdObj, std::shared_ptr<CollidingObject> cutObj);
+    virtual ~PbdObjectCuttingPair() override = default;
+
+    void apply();
+
+protected:
+
+    ///
+    /// \brief Add new vertices to pbdObj
+    ///
+    void addVertices(std::shared_ptr<SurfaceMesh> pbdMesh,
+                     std::shared_ptr<VecDataArray<double, 3>> vertices,
+                     std::shared_ptr<VecDataArray<double, 3>> initialVertices);
+
+    ///
+    /// \brief Modify current vertices of pbdObj
+    ///
+    void modifyVertices(std::shared_ptr<SurfaceMesh> pbdMesh,
+                        std::shared_ptr<std::vector<size_t>> vertexIndices,
+                        std::shared_ptr<VecDataArray<double, 3>> vertices,
+                        std::shared_ptr<VecDataArray<double, 3>> initialVertices);
+
+    ///
+    /// \brief Add new elements to pbdObj
+    ///
+    void addTriangles(std::shared_ptr<SurfaceMesh> pbdMesh,
+                      std::shared_ptr<VecDataArray<int, 3>> elements);
+
+    ///
+    /// \brief Modify exsiting elements of pbdObj
+    ///
+    void modifyTriangles(std::shared_ptr<SurfaceMesh> pbdMesh,
+                         std::shared_ptr<std::vector<size_t>> elementIndices,
+                         std::shared_ptr<VecDataArray<int, 3>> elements);
+
+    std::shared_ptr<std::unordered_set<size_t>> m_removeConstraintVertices = std::make_shared<std::unordered_set<size_t>>();
+    std::shared_ptr<std::unordered_set<size_t>> m_addConstraintVertices    = std::make_shared<std::unordered_set<size_t>>();
+};
+}
\ No newline at end of file
diff --git a/Source/Scene/imstkPbdObjectPickingPair.h b/Source/Scene/imstkPbdObjectPickingPair.h
index 6ba0ca6d34e22363520d48a3fd9a347fcf01b6f2..e5a6fc923f72d14e54b687882b98df634ca9191f 100644
--- a/Source/Scene/imstkPbdObjectPickingPair.h
+++ b/Source/Scene/imstkPbdObjectPickingPair.h
@@ -32,7 +32,7 @@ class PbdSolver;
 ///
 /// \class PbdObjectCollisionPair
 ///
-/// \brief This class defines a collision interaction between two PbdObjects
+/// \brief This class defines a picking interaction between a PbdObject and a CollidingObject with AnalyticalGeometry
 ///
 class PbdObjectPickingPair : public CollisionPair
 {
diff --git a/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.cpp b/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.cpp
index fc3bebbc147a0fc3960e1e6798eb9bcb522d34d7..7ebd5b2ad2f37946164e13fa035eb278f90f2ce7 100644
--- a/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.cpp
+++ b/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.cpp
@@ -43,16 +43,17 @@
 
 vtkStandardNewMacro(vtkXRenderWindowInteractor2);
 
-template <int EventType>
-int XEventTypeEquals(Display*, XEvent* event, XPointer)
+template<int EventType>
+int
+XEventTypeEquals(Display*, XEvent* event, XPointer)
 {
-  return event->type == EventType;
+    return event->type == EventType;
 }
 
 struct vtkXRenderWindowInteractor2Timer
 {
-  unsigned long duration;
-  timeval lastFire;
+    unsigned long duration;
+    timeval lastFire;
 };
 
 // Map between the X native id to our own integer count id.  Note this
@@ -62,79 +63,79 @@ struct vtkXRenderWindowInteractor2Timer
 class vtkXRenderWindowInteractor2Internals
 {
 public:
-  vtkXRenderWindowInteractor2Internals() { this->TimerIdCount = 1; }
-  ~vtkXRenderWindowInteractor2Internals() = default;
-
-  // duration is in milliseconds
-  int CreateLocalTimer(unsigned long duration)
-  {
-    int id = this->TimerIdCount++;
-    this->LocalToTimer[id].duration = duration;
-    gettimeofday(&this->LocalToTimer[id].lastFire, nullptr);
-    return id;
-  }
-
-  void DestroyLocalTimer(int id) { this->LocalToTimer.erase(id); }
-
-  void GetTimeToNextTimer(timeval& tv)
-  {
-    uint64_t lowestDelta = 1000000;
-    if (!this->LocalToTimer.empty())
+    vtkXRenderWindowInteractor2Internals() { this->TimerIdCount = 1; }
+    ~vtkXRenderWindowInteractor2Internals() = default;
+
+    // duration is in milliseconds
+    int CreateLocalTimer(unsigned long duration)
+    {
+        int id = this->TimerIdCount++;
+        this->LocalToTimer[id].duration = duration;
+        gettimeofday(&this->LocalToTimer[id].lastFire, nullptr);
+        return id;
+    }
+
+    void DestroyLocalTimer(int id) { this->LocalToTimer.erase(id); }
+
+    void GetTimeToNextTimer(timeval& tv)
     {
-      timeval ctv;
-      gettimeofday(&ctv, nullptr);
-      for (auto& timer : this->LocalToTimer)
-      {
-        uint64_t delta = (ctv.tv_sec - timer.second.lastFire.tv_sec) * 1000000 + ctv.tv_usec -
-          timer.second.lastFire.tv_usec;
-        if (delta < lowestDelta)
+        uint64_t lowestDelta = 1000000;
+        if (!this->LocalToTimer.empty())
         {
-          lowestDelta = delta;
+            timeval ctv;
+            gettimeofday(&ctv, nullptr);
+            for (auto& timer : this->LocalToTimer)
+            {
+                uint64_t delta = (ctv.tv_sec - timer.second.lastFire.tv_sec) * 1000000 + ctv.tv_usec -
+                                 timer.second.lastFire.tv_usec;
+                if (delta < lowestDelta)
+                {
+                    lowestDelta = delta;
+                }
+            }
         }
-      }
+        tv.tv_sec  = lowestDelta / 1000000;
+        tv.tv_usec = lowestDelta % 1000000;
     }
-    tv.tv_sec = lowestDelta / 1000000;
-    tv.tv_usec = lowestDelta % 1000000;
-  }
 
-  void FireTimers(vtkXRenderWindowInteractor2* rwi)
-  {
-    if (!this->LocalToTimer.empty())
+    void FireTimers(vtkXRenderWindowInteractor2* rwi)
     {
-      timeval ctv;
-      gettimeofday(&ctv, nullptr);
-      std::vector<unsigned long> expired;
-      for (auto& timer : this->LocalToTimer)
-      {
-        int64_t delta = (ctv.tv_sec - timer.second.lastFire.tv_sec) * 1000000 + ctv.tv_usec -
-          timer.second.lastFire.tv_usec;
-        if (delta / 1000 >= static_cast<int64_t>(timer.second.duration))
+        if (!this->LocalToTimer.empty())
         {
-          int timerId = rwi->GetVTKTimerId(timer.first);
-          rwi->InvokeEvent(vtkCommand::TimerEvent, &timerId);
-          if (rwi->IsOneShotTimer(timerId))
-          {
-            expired.push_back(timer.first);
-          }
-          else
-          {
-            timer.second.lastFire.tv_sec = ctv.tv_sec;
-            timer.second.lastFire.tv_usec = ctv.tv_usec;
-          }
+            timeval ctv;
+            gettimeofday(&ctv, nullptr);
+            std::vector<unsigned long> expired;
+            for (auto& timer : this->LocalToTimer)
+            {
+                int64_t delta = (ctv.tv_sec - timer.second.lastFire.tv_sec) * 1000000 + ctv.tv_usec -
+                                timer.second.lastFire.tv_usec;
+                if (delta / 1000 >= static_cast<int64_t>(timer.second.duration))
+                {
+                    int timerId = rwi->GetVTKTimerId(timer.first);
+                    rwi->InvokeEvent(vtkCommand::TimerEvent, &timerId);
+                    if (rwi->IsOneShotTimer(timerId))
+                    {
+                        expired.push_back(timer.first);
+                    }
+                    else
+                    {
+                        timer.second.lastFire.tv_sec  = ctv.tv_sec;
+                        timer.second.lastFire.tv_usec = ctv.tv_usec;
+                    }
+                }
+            }
+            for (auto exp : expired)
+            {
+                this->DestroyLocalTimer(exp);
+            }
         }
-      }
-      for (auto exp : expired)
-      {
-        this->DestroyLocalTimer(exp);
-      }
     }
-  }
 
-  static std::set<vtkXRenderWindowInteractor2*> Instances;
+    static std::set<vtkXRenderWindowInteractor2*> Instances;
 
 private:
-  int TimerIdCount;
-  std::map<int, vtkXRenderWindowInteractor2Timer> LocalToTimer;
+    int TimerIdCount;
+    std::map<int, vtkXRenderWindowInteractor2Timer> LocalToTimer;
 };
 
 std::set<vtkXRenderWindowInteractor2*> vtkXRenderWindowInteractor2Internals::Instances;
@@ -145,716 +146,729 @@ typedef XID vtkKeySym;
 //------------------------------------------------------------------------------
 vtkXRenderWindowInteractor2::vtkXRenderWindowInteractor2()
 {
-  this->Internal = new vtkXRenderWindowInteractor2Internals;
-  this->DisplayId = nullptr;
-  this->WindowId = 0;
-  this->KillAtom = 0;
-  this->XdndSource = 0;
-  this->XdndPositionAtom = 0;
-  this->XdndDropAtom = 0;
-  this->XdndActionCopyAtom = 0;
-  this->XdndStatusAtom = 0;
-  this->XdndFinishedAtom = 0;
+    this->Internal           = new vtkXRenderWindowInteractor2Internals;
+    this->DisplayId          = nullptr;
+    this->WindowId           = 0;
+    this->KillAtom           = 0;
+    this->XdndSource         = 0;
+    this->XdndPositionAtom   = 0;
+    this->XdndDropAtom       = 0;
+    this->XdndActionCopyAtom = 0;
+    this->XdndStatusAtom     = 0;
+    this->XdndFinishedAtom   = 0;
 }
 
 //------------------------------------------------------------------------------
 vtkXRenderWindowInteractor2::~vtkXRenderWindowInteractor2()
 {
-  this->Disable();
+    this->Disable();
 
-  delete this->Internal;
+    delete this->Internal;
 }
 
 //------------------------------------------------------------------------------
 // TerminateApp() notifies the event loop to exit.
 // The event loop is started by Start() or by one own's method.
 // This results in Start() returning to its caller.
-void vtkXRenderWindowInteractor2::TerminateApp()
+void
+vtkXRenderWindowInteractor2::TerminateApp()
 {
-  if (this->Done)
-  {
-    return;
-  }
-
-  this->Done = true;
-
-  // Send a VTK_BreakXtLoop ClientMessage event to be sure we pop out of the
-  // event loop.  This "wakes up" the event loop.  Otherwise, it might sit idle
-  // waiting for an event before realizing an exit was requested.
-  XClientMessageEvent client;
-  memset(&client, 0, sizeof(client));
-
-  client.type = ClientMessage;
-  // client.serial; //leave zeroed
-  // client.send_event; //leave zeroed
-  client.display = this->DisplayId;
-  client.window = this->WindowId;
-  client.message_type = XInternAtom(this->DisplayId, "VTK_BreakXtLoop", False);
-  client.format = 32; // indicates size of data chunks: 8, 16 or 32 bits...
-  // client.data; //leave zeroed
-
-  XSendEvent(client.display, client.window, True, NoEventMask, reinterpret_cast<XEvent*>(&client));
-  XFlush(client.display);
-  this->RenderWindow->Finalize();
+    if (this->Done)
+    {
+        return;
+    }
+
+    this->Done = true;
+
+    // Send a VTK_BreakXtLoop ClientMessage event to be sure we pop out of the
+    // event loop.  This "wakes up" the event loop.  Otherwise, it might sit idle
+    // waiting for an event before realizing an exit was requested.
+    XClientMessageEvent client;
+    memset(&client, 0, sizeof(client));
+
+    client.type = ClientMessage;
+    // client.serial; //leave zeroed
+    // client.send_event; //leave zeroed
+    client.display      = this->DisplayId;
+    client.window       = this->WindowId;
+    client.message_type = XInternAtom(this->DisplayId, "VTK_BreakXtLoop", False);
+    client.format       = 32; // indicates size of data chunks: 8, 16 or 32 bits...
+    // client.data; //leave zeroed
+
+    XSendEvent(client.display, client.window, True, NoEventMask, reinterpret_cast<XEvent*>(&client));
+    XFlush(client.display);
+    this->RenderWindow->Finalize();
 }
 
-void vtkXRenderWindowInteractor2::ProcessEvents()
+void
+vtkXRenderWindowInteractor2::ProcessEvents()
 {
-  XEvent event;
-  while (XPending(this->DisplayId) && !this->Done)
-  {
-    XNextEvent(this->DisplayId, &event);
-    this->DispatchEvent(&event);
-  }
+    XEvent event;
+    while (XPending(this->DisplayId) && !this->Done)
+    {
+        XNextEvent(this->DisplayId, &event);
+        this->DispatchEvent(&event);
+    }
 }
 
 //------------------------------------------------------------------------------
 // This will start up the X event loop. If you
 // call this method it will loop processing X events until the
 // loop is exited.
-void vtkXRenderWindowInteractor2::StartEventLoop()
+void
+vtkXRenderWindowInteractor2::StartEventLoop()
 {
-  std::vector<int> rwiFileDescriptors;
-  fd_set in_fds;
-  struct timeval tv;
-  struct timeval minTv;
-
-  for (auto rwi : vtkXRenderWindowInteractor2Internals::Instances)
-  {
-    rwi->Done = false;
-    rwiFileDescriptors.push_back(ConnectionNumber(rwi->DisplayId));
-  }
-
-  bool done = true;
-  do
-  {
-    bool wait = true;
-    done = true;
-    minTv.tv_sec = 1000;
-    minTv.tv_usec = 1000;
-    XEvent event;
-    for (auto rwi = vtkXRenderWindowInteractor2Internals::Instances.begin();
-         rwi != vtkXRenderWindowInteractor2Internals::Instances.end();)
+    std::vector<int> rwiFileDescriptors;
+    fd_set           in_fds;
+    struct timeval   tv;
+    struct timeval   minTv;
+
+    for (auto rwi : vtkXRenderWindowInteractor2Internals::Instances)
     {
+        rwi->Done = false;
+        rwiFileDescriptors.push_back(ConnectionNumber(rwi->DisplayId));
+    }
 
-      if (XPending((*rwi)->DisplayId) == 0)
-      {
-        // get how long to wait for the next timer
-        (*rwi)->Internal->GetTimeToNextTimer(tv);
-        minTv.tv_sec = std::min(tv.tv_sec, minTv.tv_sec);
-        minTv.tv_usec = std::min(tv.tv_usec, minTv.tv_usec);
-      }
-      else
-      {
-        // If events are pending, dispatch them to the right RenderWindowInteractor
-        XNextEvent((*rwi)->DisplayId, &event);
-        (*rwi)->DispatchEvent(&event);
-        wait = false;
-      }
-      (*rwi)->FireTimers();
-
-      // Check if all RenderWindowInteractors have been terminated
-      done = done && (*rwi)->Done;
-
-      // If current RenderWindowInteractor have been terminated, handle its last event,
-      // then remove it from the Instance vector
-      if ((*rwi)->Done)
-      {
-        // Empty the event list
-        while (XPending((*rwi)->DisplayId) != 0)
+    bool done = true;
+    do
+    {
+        bool wait = true;
+        done = true;
+        minTv.tv_sec  = 1000;
+        minTv.tv_usec = 1000;
+        XEvent event;
+        for (auto rwi = vtkXRenderWindowInteractor2Internals::Instances.begin();
+             rwi != vtkXRenderWindowInteractor2Internals::Instances.end();)
         {
-          XNextEvent((*rwi)->DisplayId, &event);
-          (*rwi)->DispatchEvent(&event);
+            if (XPending((*rwi)->DisplayId) == 0)
+            {
+                // get how long to wait for the next timer
+                (*rwi)->Internal->GetTimeToNextTimer(tv);
+                minTv.tv_sec  = std::min(tv.tv_sec, minTv.tv_sec);
+                minTv.tv_usec = std::min(tv.tv_usec, minTv.tv_usec);
+            }
+            else
+            {
+                // If events are pending, dispatch them to the right RenderWindowInteractor
+                XNextEvent((*rwi)->DisplayId, &event);
+                (*rwi)->DispatchEvent(&event);
+                wait = false;
+            }
+            (*rwi)->FireTimers();
+
+            // Check if all RenderWindowInteractors have been terminated
+            done = done && (*rwi)->Done;
+
+            // If current RenderWindowInteractor have been terminated, handle its last event,
+            // then remove it from the Instance vector
+            if ((*rwi)->Done)
+            {
+                // Empty the event list
+                while (XPending((*rwi)->DisplayId) != 0)
+                {
+                    XNextEvent((*rwi)->DisplayId, &event);
+                    (*rwi)->DispatchEvent(&event);
+                }
+                // Adjust the file descriptors vector
+                int rwiPosition =
+                    std::distance(vtkXRenderWindowInteractor2Internals::Instances.begin(), rwi);
+                rwi = vtkXRenderWindowInteractor2Internals::Instances.erase(rwi);
+                rwiFileDescriptors.erase(rwiFileDescriptors.begin() + rwiPosition);
+            }
+            else
+            {
+                ++rwi;
+            }
         }
-        // Adjust the file descriptors vector
-        int rwiPosition =
-          std::distance(vtkXRenderWindowInteractor2Internals::Instances.begin(), rwi);
-        rwi = vtkXRenderWindowInteractor2Internals::Instances.erase(rwi);
-        rwiFileDescriptors.erase(rwiFileDescriptors.begin() + rwiPosition);
-      }
-      else
-      {
-        ++rwi;
-      }
-    }
 
-    if (wait && !done)
-    {
-      // select will wait until 'tv' elapses or something else wakes us
-      FD_ZERO(&in_fds);
-      for (auto rwiFileDescriptor : rwiFileDescriptors)
-      {
-        FD_SET(rwiFileDescriptor, &in_fds);
-      }
-      int maxFileDescriptor =
-        *std::max_element(rwiFileDescriptors.begin(), rwiFileDescriptors.end());
-      select(maxFileDescriptor + 1, &in_fds, nullptr, nullptr, &minTv);
+        if (wait && !done)
+        {
+            // select will wait until 'tv' elapses or something else wakes us
+            FD_ZERO(&in_fds);
+            for (auto rwiFileDescriptor : rwiFileDescriptors)
+            {
+                FD_SET(rwiFileDescriptor, &in_fds);
+            }
+            int maxFileDescriptor =
+                *std::max_element(rwiFileDescriptors.begin(), rwiFileDescriptors.end());
+            select(maxFileDescriptor + 1, &in_fds, nullptr, nullptr, &minTv);
+        }
     }
-
-  } while (!done);
+    while (!done);
 }
 
 //------------------------------------------------------------------------------
 // Initializes the event handlers without an XtAppContext.  This is
 // good for when you don't have a user interface, but you still
 // want to have mouse interaction.
-void vtkXRenderWindowInteractor2::Initialize()
+void
+vtkXRenderWindowInteractor2::Initialize()
 {
-  if (this->Initialized)
-  {
-    return;
-  }
-
-  vtkRenderWindow* ren;
-  int* size;
-
-  // make sure we have a RenderWindow and camera
-  if (!this->RenderWindow)
-  {
-    vtkErrorMacro(<< "No renderer defined!");
-    return;
-  }
-
-  this->Initialized = 1;
-  ren = this->RenderWindow;
-
-  this->DisplayId = static_cast<Display*>(ren->GetGenericDisplayId());
-  if (!this->DisplayId)
-  {
-    vtkDebugMacro("opening display");
-    this->DisplayId = XOpenDisplay(nullptr);
-    vtkDebugMacro("opened display");
-    ren->SetDisplayId(this->DisplayId);
-  }
-
-  vtkXRenderWindowInteractor2Internals::Instances.insert(this);
-
-  size = ren->GetActualSize();
-  size[0] = ((size[0] > 0) ? size[0] : 300);
-  size[1] = ((size[1] > 0) ? size[1] : 300);
-  XSync(this->DisplayId, False);
-
-  ren->Start();
-  ren->End();
-
-  this->WindowId = reinterpret_cast<Window>(ren->GetGenericWindowId());
-
-  XWindowAttributes attribs;
-  //  Find the current window size
-  XGetWindowAttributes(this->DisplayId, this->WindowId, &attribs);
-
-  size[0] = attribs.width;
-  size[1] = attribs.height;
-  ren->SetSize(size[0], size[1]);
-
-  this->Enable();
-  this->Size[0] = size[0];
-  this->Size[1] = size[1];
+    if (this->Initialized)
+    {
+        return;
+    }
+
+    vtkRenderWindow* ren;
+    int*             size;
+
+    // make sure we have a RenderWindow and camera
+    if (!this->RenderWindow)
+    {
+        vtkErrorMacro(<< "No renderer defined!");
+        return;
+    }
+
+    this->Initialized = 1;
+    ren = this->RenderWindow;
+
+    this->DisplayId = static_cast<Display*>(ren->GetGenericDisplayId());
+    if (!this->DisplayId)
+    {
+        vtkDebugMacro("opening display");
+        this->DisplayId = XOpenDisplay(nullptr);
+        vtkDebugMacro("opened display");
+        ren->SetDisplayId(this->DisplayId);
+    }
+
+    vtkXRenderWindowInteractor2Internals::Instances.insert(this);
+
+    size    = ren->GetActualSize();
+    size[0] = ((size[0] > 0) ? size[0] : 300);
+    size[1] = ((size[1] > 0) ? size[1] : 300);
+    XSync(this->DisplayId, False);
+
+    ren->Start();
+    ren->End();
+
+    this->WindowId = reinterpret_cast<Window>(ren->GetGenericWindowId());
+
+    XWindowAttributes attribs;
+    //  Find the current window size
+    XGetWindowAttributes(this->DisplayId, this->WindowId, &attribs);
+
+    size[0] = attribs.width;
+    size[1] = attribs.height;
+    ren->SetSize(size[0], size[1]);
+
+    this->Enable();
+    this->Size[0] = size[0];
+    this->Size[1] = size[1];
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::Enable()
+void
+vtkXRenderWindowInteractor2::Enable()
 {
-  // avoid cycles of calling Initialize() and Enable()
-  if (this->Enabled)
-  {
-    return;
-  }
-
-  // Add the event handler to the system.
-  // If we change the types of events processed by this handler, then
-  // we need to change the Disable() routine to match.  In order for Disable()
-  // to work properly, both the callback function AND the client data
-  // passed to XtAddEventHandler and XtRemoveEventHandler must MATCH
-  // PERFECTLY
-  XSelectInput(this->DisplayId, this->WindowId,
+    // avoid cycles of calling Initialize() and Enable()
+    if (this->Enabled)
+    {
+        return;
+    }
+
+    // Add the event handler to the system.
+    // If we change the types of events processed by this handler, then
+    // we need to change the Disable() routine to match.  In order for Disable()
+    // to work properly, both the callback function AND the client data
+    // passed to XtAddEventHandler and XtRemoveEventHandler must MATCH
+    // PERFECTLY
+    XSelectInput(this->DisplayId, this->WindowId,
     KeyPressMask | KeyReleaseMask | ButtonPressMask | ButtonReleaseMask | ExposureMask |
       StructureNotifyMask | EnterWindowMask | LeaveWindowMask | PointerMotionHintMask |
       PointerMotionMask);
 
-  // Setup for capturing the window deletion
-  this->KillAtom = XInternAtom(this->DisplayId, "WM_DELETE_WINDOW", False);
-  XSetWMProtocols(this->DisplayId, this->WindowId, &this->KillAtom, 1);
+    // Setup for capturing the window deletion
+    this->KillAtom = XInternAtom(this->DisplayId, "WM_DELETE_WINDOW", False);
+    XSetWMProtocols(this->DisplayId, this->WindowId, &this->KillAtom, 1);
 
-  // Enable drag and drop
-  Atom xdndAwareAtom = XInternAtom(this->DisplayId, "XdndAware", False);
-  char xdndVersion = 5;
-  XChangeProperty(this->DisplayId, this->WindowId, xdndAwareAtom, XA_ATOM, 32, PropModeReplace,
-    (unsigned char*)&xdndVersion, 1);
-  this->XdndPositionAtom = XInternAtom(this->DisplayId, "XdndPosition", False);
-  this->XdndDropAtom = XInternAtom(this->DisplayId, "XdndDrop", False);
-  this->XdndActionCopyAtom = XInternAtom(this->DisplayId, "XdndActionCopy", False);
-  this->XdndStatusAtom = XInternAtom(this->DisplayId, "XdndStatus", False);
-  this->XdndFinishedAtom = XInternAtom(this->DisplayId, "XdndFinished", False);
+    // Enable drag and drop
+    Atom xdndAwareAtom = XInternAtom(this->DisplayId, "XdndAware", False);
+    char xdndVersion   = 5;
+    XChangeProperty(this->DisplayId, this->WindowId, xdndAwareAtom, XA_ATOM, 32, PropModeReplace,
+        (unsigned char*)&xdndVersion, 1);
+    this->XdndPositionAtom   = XInternAtom(this->DisplayId, "XdndPosition", False);
+    this->XdndDropAtom       = XInternAtom(this->DisplayId, "XdndDrop", False);
+    this->XdndActionCopyAtom = XInternAtom(this->DisplayId, "XdndActionCopy", False);
+    this->XdndStatusAtom     = XInternAtom(this->DisplayId, "XdndStatus", False);
+    this->XdndFinishedAtom   = XInternAtom(this->DisplayId, "XdndFinished", False);
 
-  this->Enabled = 1;
+    this->Enabled = 1;
 
-  this->Modified();
+    this->Modified();
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::Disable()
+void
+vtkXRenderWindowInteractor2::Disable()
 {
-  if (!this->Enabled)
-  {
-    return;
-  }
+    if (!this->Enabled)
+    {
+        return;
+    }
 
-  this->Enabled = 0;
+    this->Enabled = 0;
 
-  this->Modified();
+    this->Modified();
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::PrintSelf(ostream& os, vtkIndent indent)
+void
+vtkXRenderWindowInteractor2::PrintSelf(ostream& os, vtkIndent indent)
 {
-  this->Superclass::PrintSelf(os, indent);
+    this->Superclass::PrintSelf(os, indent);
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::UpdateSize(int x, int y)
+void
+vtkXRenderWindowInteractor2::UpdateSize(int x, int y)
 {
-  // if the size changed send this on to the RenderWindow
-  if ((x != this->Size[0]) || (y != this->Size[1]))
-  {
-    this->Size[0] = x;
-    this->Size[1] = y;
-    this->RenderWindow->SetSize(x, y);
-  }
+    // if the size changed send this on to the RenderWindow
+    if ((x != this->Size[0]) || (y != this->Size[1]))
+    {
+        this->Size[0] = x;
+        this->Size[1] = y;
+        this->RenderWindow->SetSize(x, y);
+    }
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::UpdateSizeNoXResize(int x, int y)
+void
+vtkXRenderWindowInteractor2::UpdateSizeNoXResize(int x, int y)
 {
-  // if the size changed send this on to the RenderWindow
-  if ((x != this->Size[0]) || (y != this->Size[1]))
-  {
-    this->Size[0] = x;
-    this->Size[1] = y;
-    // static_cast<vtkXOpenGLRenderWindow*>(this->RenderWindow)->SetSizeNoXResize(x, y);
-    this->RenderWindow->SetSize(x, y);
-  }
+    // if the size changed send this on to the RenderWindow
+    if ((x != this->Size[0]) || (y != this->Size[1]))
+    {
+        this->Size[0] = x;
+        this->Size[1] = y;
+        // static_cast<vtkXOpenGLRenderWindow*>(this->RenderWindow)->SetSizeNoXResize(x, y);
+        this->RenderWindow->SetSize(x, y);
+    }
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::FireTimers()
+void
+vtkXRenderWindowInteractor2::FireTimers()
 {
-  if (this->GetEnabled())
-  {
-    this->Internal->FireTimers(this);
-  }
+    if (this->GetEnabled())
+    {
+        this->Internal->FireTimers(this);
+    }
 }
 
 //------------------------------------------------------------------------------
 // X always creates one shot timers
-int vtkXRenderWindowInteractor2::InternalCreateTimer(
-  int vtkNotUsed(timerId), int vtkNotUsed(timerType), unsigned long duration)
+int
+vtkXRenderWindowInteractor2::InternalCreateTimer(
+    int vtkNotUsed(timerId), int vtkNotUsed(timerType), unsigned long duration)
 {
-  duration = (duration > 0 ? duration : this->TimerDuration);
-  int platformTimerId = this->Internal->CreateLocalTimer(duration);
-  return platformTimerId;
+    duration = (duration > 0 ? duration : this->TimerDuration);
+    int platformTimerId = this->Internal->CreateLocalTimer(duration);
+    return platformTimerId;
 }
 
 //------------------------------------------------------------------------------
-int vtkXRenderWindowInteractor2::InternalDestroyTimer(int platformTimerId)
+int
+vtkXRenderWindowInteractor2::InternalDestroyTimer(int platformTimerId)
 {
-  this->Internal->DestroyLocalTimer(platformTimerId);
-  return 1;
+    this->Internal->DestroyLocalTimer(platformTimerId);
+    return 1;
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::DispatchEvent(XEvent* event)
+void
+vtkXRenderWindowInteractor2::DispatchEvent(XEvent* event)
 {
-  int xp, yp;
+    int xp, yp;
 
-  switch (event->type)
-  {
+    switch (event->type)
+    {
     case Expose:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      XEvent result;
-      while (XCheckTypedWindowEvent(this->DisplayId, this->WindowId, Expose, &result))
-      {
-        // just getting the expose configure event
-        event = &result;
-      }
-      XExposeEvent* exposeEvent = reinterpret_cast<XExposeEvent*>(event);
-      this->SetEventSize(exposeEvent->width, exposeEvent->height);
-      xp = exposeEvent->x;
-      yp = exposeEvent->y;
-      yp = this->Size[1] - yp - 1;
-      this->SetEventPosition(xp, yp);
-
-      // only render if we are currently accepting events
-      if (this->Enabled)
-      {
-        this->InvokeEvent(vtkCommand::ExposeEvent, nullptr);
-        this->Render();
-      }
+        if (!this->Enabled)
+        {
+            return;
+        }
+        XEvent result;
+        while (XCheckTypedWindowEvent(this->DisplayId, this->WindowId, Expose, &result))
+        {
+            // just getting the expose configure event
+            event = &result;
+        }
+        XExposeEvent* exposeEvent = reinterpret_cast<XExposeEvent*>(event);
+        this->SetEventSize(exposeEvent->width, exposeEvent->height);
+        xp = exposeEvent->x;
+        yp = exposeEvent->y;
+        yp = this->Size[1] - yp - 1;
+        this->SetEventPosition(xp, yp);
+
+        // only render if we are currently accepting events
+        if (this->Enabled)
+        {
+            this->InvokeEvent(vtkCommand::ExposeEvent, nullptr);
+            this->Render();
+        }
     }
     break;
 
     case MapNotify:
     {
-      // only render if we are currently accepting events
-      if (this->Enabled && this->GetRenderWindow()->GetNeverRendered())
-      {
-        this->Render();
-      }
+        // only render if we are currently accepting events
+        if (this->Enabled && this->GetRenderWindow()->GetNeverRendered())
+        {
+            this->Render();
+        }
     }
     break;
 
     case ConfigureNotify:
     {
-      XEvent result;
-      while (XCheckTypedWindowEvent(this->DisplayId, this->WindowId, ConfigureNotify, &result))
-      {
-        // just getting the last configure event
-        event = &result;
-      }
-      int width = (reinterpret_cast<XConfigureEvent*>(event))->width;
-      int height = (reinterpret_cast<XConfigureEvent*>(event))->height;
-      if (width != this->Size[0] || height != this->Size[1])
-      {
-        bool resizeSmaller = width <= this->Size[0] && height <= this->Size[1];
-        this->UpdateSizeNoXResize(width, height);
-        xp = (reinterpret_cast<XButtonEvent*>(event))->x;
-        yp = (reinterpret_cast<XButtonEvent*>(event))->y;
-        this->SetEventPosition(xp, this->Size[1] - yp - 1);
-        // only render if we are currently accepting events
-        if (this->Enabled)
+        XEvent result;
+        while (XCheckTypedWindowEvent(this->DisplayId, this->WindowId, ConfigureNotify, &result))
         {
-          this->InvokeEvent(vtkCommand::ConfigureEvent, nullptr);
-          if (resizeSmaller)
-          {
-            // Don't call Render when the window is resized to be larger:
-            //
-            // - if the window is resized to be larger, an Expose event will
-            // be triggered by the X server which will trigger a call to
-            // Render().
-            // - if the window is resized to be smaller, no Expose event will
-            // be triggered by the X server, as no new area become visible.
-            // only in this case, we need to explicitly call Render()
-            // in ConfigureNotify.
-            this->Render();
-          }
+            // just getting the last configure event
+            event = &result;
+        }
+        int width  = (reinterpret_cast<XConfigureEvent*>(event))->width;
+        int height = (reinterpret_cast<XConfigureEvent*>(event))->height;
+        if (width != this->Size[0] || height != this->Size[1])
+        {
+            bool resizeSmaller = width <= this->Size[0] && height <= this->Size[1];
+            this->UpdateSizeNoXResize(width, height);
+            xp = (reinterpret_cast<XButtonEvent*>(event))->x;
+            yp = (reinterpret_cast<XButtonEvent*>(event))->y;
+            this->SetEventPosition(xp, this->Size[1] - yp - 1);
+            // only render if we are currently accepting events
+            if (this->Enabled)
+            {
+                this->InvokeEvent(vtkCommand::ConfigureEvent, nullptr);
+                if (resizeSmaller)
+                {
+                    // Don't call Render when the window is resized to be larger:
+                    //
+                    // - if the window is resized to be larger, an Expose event will
+                    // be triggered by the X server which will trigger a call to
+                    // Render().
+                    // - if the window is resized to be smaller, no Expose event will
+                    // be triggered by the X server, as no new area become visible.
+                    // only in this case, we need to explicitly call Render()
+                    // in ConfigureNotify.
+                    this->Render();
+                }
+            }
         }
-      }
     }
     break;
 
     case ButtonPress:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      int ctrl = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
-      int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
-      int alt = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
-      xp = (reinterpret_cast<XButtonEvent*>(event))->x;
-      yp = (reinterpret_cast<XButtonEvent*>(event))->y;
-
-      // check for double click
-      static int MousePressTime = 0;
-      int repeat = 0;
-      // 400 ms threshold by default is probably good to start
-      int eventTime = static_cast<int>(reinterpret_cast<XButtonEvent*>(event)->time);
-      if ((eventTime - MousePressTime) < 400)
-      {
-        MousePressTime -= 2000; // no double click next time
-        repeat = 1;
-      }
-      else
-      {
-        MousePressTime = eventTime;
-      }
-
-      this->SetEventInformationFlipY(xp, yp, ctrl, shift, 0, repeat);
-      this->SetAltKey(alt);
-      switch ((reinterpret_cast<XButtonEvent*>(event))->button)
-      {
+        if (!this->Enabled)
+        {
+            return;
+        }
+        int ctrl  = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
+        int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
+        int alt   = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
+        xp = (reinterpret_cast<XButtonEvent*>(event))->x;
+        yp = (reinterpret_cast<XButtonEvent*>(event))->y;
+
+        // check for double click
+        static int MousePressTime = 0;
+        int        repeat = 0;
+        // 400 ms threshold by default is probably good to start
+        int eventTime = static_cast<int>(reinterpret_cast<XButtonEvent*>(event)->time);
+        if ((eventTime - MousePressTime) < 400)
+        {
+            MousePressTime -= 2000; // no double click next time
+            repeat = 1;
+        }
+        else
+        {
+            MousePressTime = eventTime;
+        }
+
+        this->SetEventInformationFlipY(xp, yp, ctrl, shift, 0, repeat);
+        this->SetAltKey(alt);
+        switch ((reinterpret_cast<XButtonEvent*>(event))->button)
+        {
         case Button1:
-          this->InvokeEvent(vtkCommand::LeftButtonPressEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::LeftButtonPressEvent, nullptr);
+            break;
         case Button2:
-          this->InvokeEvent(vtkCommand::MiddleButtonPressEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::MiddleButtonPressEvent, nullptr);
+            break;
         case Button3:
-          this->InvokeEvent(vtkCommand::RightButtonPressEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::RightButtonPressEvent, nullptr);
+            break;
         case Button4:
-          this->InvokeEvent(vtkCommand::MouseWheelForwardEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::MouseWheelForwardEvent, nullptr);
+            break;
         case Button5:
-          this->InvokeEvent(vtkCommand::MouseWheelBackwardEvent, nullptr);
-          break;
-      }
+            this->InvokeEvent(vtkCommand::MouseWheelBackwardEvent, nullptr);
+            break;
+        }
     }
     break;
 
     case ButtonRelease:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      int ctrl = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
-      int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
-      int alt = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
-      xp = (reinterpret_cast<XButtonEvent*>(event))->x;
-      yp = (reinterpret_cast<XButtonEvent*>(event))->y;
-      this->SetEventInformationFlipY(xp, yp, ctrl, shift);
-      this->SetAltKey(alt);
-      switch ((reinterpret_cast<XButtonEvent*>(event))->button)
-      {
+        if (!this->Enabled)
+        {
+            return;
+        }
+        int ctrl  = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
+        int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
+        int alt   = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
+        xp = (reinterpret_cast<XButtonEvent*>(event))->x;
+        yp = (reinterpret_cast<XButtonEvent*>(event))->y;
+        this->SetEventInformationFlipY(xp, yp, ctrl, shift);
+        this->SetAltKey(alt);
+        switch ((reinterpret_cast<XButtonEvent*>(event))->button)
+        {
         case Button1:
-          this->InvokeEvent(vtkCommand::LeftButtonReleaseEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::LeftButtonReleaseEvent, nullptr);
+            break;
         case Button2:
-          this->InvokeEvent(vtkCommand::MiddleButtonReleaseEvent, nullptr);
-          break;
+            this->InvokeEvent(vtkCommand::MiddleButtonReleaseEvent, nullptr);
+            break;
         case Button3:
-          this->InvokeEvent(vtkCommand::RightButtonReleaseEvent, nullptr);
-          break;
-      }
+            this->InvokeEvent(vtkCommand::RightButtonReleaseEvent, nullptr);
+            break;
+        }
     }
     break;
 
     case EnterNotify:
     {
-      // Force the keyboard focus to be this render window
-      XSetInputFocus(this->DisplayId, this->WindowId, RevertToPointerRoot, CurrentTime);
-      if (this->Enabled)
-      {
-        XEnterWindowEvent* e = reinterpret_cast<XEnterWindowEvent*>(event);
-        this->SetEventInformationFlipY(
+        // Force the keyboard focus to be this render window
+        XSetInputFocus(this->DisplayId, this->WindowId, RevertToPointerRoot, CurrentTime);
+        if (this->Enabled)
+        {
+            XEnterWindowEvent* e = reinterpret_cast<XEnterWindowEvent*>(event);
+            this->SetEventInformationFlipY(
           e->x, e->y, (e->state & ControlMask) != 0, (e->state & ShiftMask) != 0);
-        this->SetAltKey(((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0);
-        this->InvokeEvent(vtkCommand::EnterEvent, nullptr);
-      }
+            this->SetAltKey(((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0);
+            this->InvokeEvent(vtkCommand::EnterEvent, nullptr);
+        }
     }
     break;
 
     case LeaveNotify:
     {
-      if (this->Enabled)
-      {
-        XLeaveWindowEvent* e = reinterpret_cast<XLeaveWindowEvent*>(event);
-        this->SetEventInformationFlipY(
+        if (this->Enabled)
+        {
+            XLeaveWindowEvent* e = reinterpret_cast<XLeaveWindowEvent*>(event);
+            this->SetEventInformationFlipY(
           e->x, e->y, (e->state & ControlMask) != 0, (e->state & ShiftMask) != 0);
-        this->SetAltKey(((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0);
-        this->InvokeEvent(vtkCommand::LeaveEvent, nullptr);
-      }
+            this->SetAltKey(((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0);
+            this->InvokeEvent(vtkCommand::LeaveEvent, nullptr);
+        }
     }
     break;
 
     case KeyPress:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      int ctrl = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
-      int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
-      int alt = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
-      vtkKeySym ks;
-      static char buffer[20];
-      buffer[0] = '\0';
-      XLookupString(reinterpret_cast<XKeyEvent*>(event), buffer, 20, &ks, nullptr);
-      xp = (reinterpret_cast<XKeyEvent*>(event))->x;
-      yp = (reinterpret_cast<XKeyEvent*>(event))->y;
-      this->SetEventInformationFlipY(xp, yp, ctrl, shift, buffer[0], 1, XKeysymToString(ks));
-      this->SetAltKey(alt);
-      this->InvokeEvent(vtkCommand::KeyPressEvent, nullptr);
-      this->InvokeEvent(vtkCommand::CharEvent, nullptr);
+        if (!this->Enabled)
+        {
+            return;
+        }
+        int         ctrl  = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
+        int         shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
+        int         alt   = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
+        vtkKeySym   ks;
+        static char buffer[20];
+        buffer[0] = '\0';
+        XLookupString(reinterpret_cast<XKeyEvent*>(event), buffer, 20, &ks, nullptr);
+        xp = (reinterpret_cast<XKeyEvent*>(event))->x;
+        yp = (reinterpret_cast<XKeyEvent*>(event))->y;
+        this->SetEventInformationFlipY(xp, yp, ctrl, shift, buffer[0], 1, XKeysymToString(ks));
+        this->SetAltKey(alt);
+        this->InvokeEvent(vtkCommand::KeyPressEvent, nullptr);
+        this->InvokeEvent(vtkCommand::CharEvent, nullptr);
     }
     break;
 
     case KeyRelease:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      int ctrl = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
-      int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
-      int alt = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
-      vtkKeySym ks;
-      static char buffer[20];
-      buffer[0] = '\0';
-      XLookupString(reinterpret_cast<XKeyEvent*>(event), buffer, 20, &ks, nullptr);
-      xp = (reinterpret_cast<XKeyEvent*>(event))->x;
-      yp = (reinterpret_cast<XKeyEvent*>(event))->y;
-      this->SetEventInformationFlipY(xp, yp, ctrl, shift, buffer[0], 1, XKeysymToString(ks));
-      this->SetAltKey(alt);
-      this->InvokeEvent(vtkCommand::KeyReleaseEvent, nullptr);
+        if (!this->Enabled)
+        {
+            return;
+        }
+        int         ctrl  = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
+        int         shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
+        int         alt   = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
+        vtkKeySym   ks;
+        static char buffer[20];
+        buffer[0] = '\0';
+        XLookupString(reinterpret_cast<XKeyEvent*>(event), buffer, 20, &ks, nullptr);
+        xp = (reinterpret_cast<XKeyEvent*>(event))->x;
+        yp = (reinterpret_cast<XKeyEvent*>(event))->y;
+        this->SetEventInformationFlipY(xp, yp, ctrl, shift, buffer[0], 1, XKeysymToString(ks));
+        this->SetAltKey(alt);
+        this->InvokeEvent(vtkCommand::KeyReleaseEvent, nullptr);
     }
     break;
 
     case MotionNotify:
     {
-      if (!this->Enabled)
-      {
-        return;
-      }
-      int ctrl = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
-      int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
-      int alt = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
-
-      // Note that even though the (x,y) location of the pointer is event structure,
-      // we must call XQueryPointer for the hints (motion event compression) to
-      // work properly.
-      this->GetMousePosition(&xp, &yp);
-      this->SetEventInformation(xp, yp, ctrl, shift);
-      this->SetAltKey(alt);
-      this->InvokeEvent(vtkCommand::MouseMoveEvent, nullptr);
+        if (!this->Enabled)
+        {
+            return;
+        }
+        int ctrl  = ((reinterpret_cast<XButtonEvent*>(event))->state & ControlMask) ? 1 : 0;
+        int shift = ((reinterpret_cast<XButtonEvent*>(event))->state & ShiftMask) ? 1 : 0;
+        int alt   = ((reinterpret_cast<XButtonEvent*>(event))->state & Mod1Mask) ? 1 : 0;
+
+        // Note that even though the (x,y) location of the pointer is event structure,
+        // we must call XQueryPointer for the hints (motion event compression) to
+        // work properly.
+        this->GetMousePosition(&xp, &yp);
+        this->SetEventInformation(xp, yp, ctrl, shift);
+        this->SetAltKey(alt);
+        this->InvokeEvent(vtkCommand::MouseMoveEvent, nullptr);
     }
     break;
 
     // Selection request for drag and drop has been delivered
     case SelectionNotify:
     {
-      // Sanity checks
-      if (!event->xselection.property || !this->XdndSource)
-      {
-        return;
-      }
-
-      // Recover the dropped file
-      char* data = nullptr;
-      Atom actualType;
-      int actualFormat;
-      unsigned long itemCount, bytesAfter;
-      XGetWindowProperty(this->DisplayId, event->xselection.requestor, event->xselection.property,
+        // Sanity checks
+        if (!event->xselection.property || !this->XdndSource)
+        {
+            return;
+        }
+
+        // Recover the dropped file
+        char*         data = nullptr;
+        Atom          actualType;
+        int           actualFormat;
+        unsigned long itemCount, bytesAfter;
+        XGetWindowProperty(this->DisplayId, event->xselection.requestor, event->xselection.property,
         0, LONG_MAX, False, event->xselection.target, &actualType, &actualFormat, &itemCount,
         &bytesAfter, (unsigned char**)&data);
 
-      // Conversion checks
-      if ((event->xselection.target != AnyPropertyType && actualType != event->xselection.target) ||
-        itemCount == 0)
-      {
-        return;
-      }
-
-      // Recover filepaths from uris and invoke DropFilesEvent
-      std::stringstream uris(data);
-      std::string uri, protocol, hostname, filePath;
-      std::string unused0, unused1, unused2, unused3;
-      vtkNew<vtkStringArray> filePaths;
-      while (std::getline(uris, uri, '\n'))
-      {
-        if (vtksys::SystemTools::ParseURL(
-              uri, protocol, unused0, unused1, hostname, unused3, filePath, true))
+        // Conversion checks
+        if ((event->xselection.target != AnyPropertyType && actualType != event->xselection.target)
+            || itemCount == 0)
+        {
+            return;
+        }
+
+        // Recover filepaths from uris and invoke DropFilesEvent
+        std::stringstream      uris(data);
+        std::string            uri, protocol, hostname, filePath;
+        std::string            unused0, unused1, unused2, unused3;
+        vtkNew<vtkStringArray> filePaths;
+        while (std::getline(uris, uri, '\n'))
         {
-          if (protocol == "file" && (hostname.empty() || hostname == "localhost"))
-          {
-            // The uris can be crlf delimited, remove ending \r if any
-            if (filePath.back() == '\r')
+            if (vtksys::SystemTools::ParseURL(
+              uri, protocol, unused0, unused1, hostname, unused3, filePath, true))
             {
-              filePath.pop_back();
+                if (protocol == "file" && (hostname.empty() || hostname == "localhost"))
+                {
+                    // The uris can be crlf delimited, remove ending \r if any
+                    if (filePath.back() == '\r')
+                    {
+                        filePath.pop_back();
+                    }
+
+                    // The extracted filepath miss the first slash
+                    filePath.insert(0, "/");
+
+                    filePaths->InsertNextValue(filePath);
+                }
             }
+        }
+        this->InvokeEvent(vtkCommand::DropFilesEvent, filePaths);
+        XFree(data);
 
-            // The extracted filepath miss the first slash
-            filePath.insert(0, "/");
+        // Inform the source the the drag and drop operation was sucessfull
+        XEvent reply;
+        memset(&reply, 0, sizeof(reply));
 
-            filePaths->InsertNextValue(filePath);
-          }
-        }
-      }
-      this->InvokeEvent(vtkCommand::DropFilesEvent, filePaths);
-      XFree(data);
-
-      // Inform the source the the drag and drop operation was sucessfull
-      XEvent reply;
-      memset(&reply, 0, sizeof(reply));
-
-      reply.type = ClientMessage;
-      reply.xclient.window = event->xclient.data.l[0];
-      reply.xclient.message_type = this->XdndFinishedAtom;
-      reply.xclient.format = 32;
-      reply.xclient.data.l[0] = this->WindowId;
-      reply.xclient.data.l[1] = itemCount;
-      reply.xclient.data.l[2] = this->XdndActionCopyAtom;
-
-      XSendEvent(this->DisplayId, this->XdndSource, False, NoEventMask, &reply);
-      XFlush(this->DisplayId);
-      this->XdndSource = 0;
+        reply.type = ClientMessage;
+        reply.xclient.window       = event->xclient.data.l[0];
+        reply.xclient.message_type = this->XdndFinishedAtom;
+        reply.xclient.format       = 32;
+        reply.xclient.data.l[0]    = this->WindowId;
+        reply.xclient.data.l[1]    = itemCount;
+        reply.xclient.data.l[2]    = this->XdndActionCopyAtom;
+
+        XSendEvent(this->DisplayId, this->XdndSource, False, NoEventMask, &reply);
+        XFlush(this->DisplayId);
+        this->XdndSource = 0;
     }
     break;
 
     case ClientMessage:
     {
-      if (event->xclient.message_type == this->XdndPositionAtom)
-      {
-        // Drag and drop event inside the window
-
-        // Recover the position
-        int xWindow, yWindow;
-        int xRoot = event->xclient.data.l[2] >> 16;
-        int yRoot = event->xclient.data.l[2] & 0xffff;
-        Window root = DefaultRootWindow(this->DisplayId);
-        Window child;
-        XTranslateCoordinates(
+        if (event->xclient.message_type == this->XdndPositionAtom)
+        {
+            // Drag and drop event inside the window
+
+            // Recover the position
+            int    xWindow, yWindow;
+            int    xRoot = event->xclient.data.l[2] >> 16;
+            int    yRoot = event->xclient.data.l[2] & 0xffff;
+            Window root  = DefaultRootWindow(this->DisplayId);
+            Window child;
+            XTranslateCoordinates(
           this->DisplayId, root, this->WindowId, xRoot, yRoot, &xWindow, &yWindow, &child);
 
-        // Convert it to VTK compatible location
-        double location[2];
-        location[0] = static_cast<double>(xWindow);
-        location[1] = static_cast<double>(this->Size[1] - yWindow - 1);
-        this->InvokeEvent(vtkCommand::UpdateDropLocationEvent, location);
-
-        // Reply that we are ready to copy the dragged data
-        XEvent reply;
-        memset(&reply, 0, sizeof(reply));
+            // Convert it to VTK compatible location
+            double location[2];
+            location[0] = static_cast<double>(xWindow);
+            location[1] = static_cast<double>(this->Size[1] - yWindow - 1);
+            this->InvokeEvent(vtkCommand::UpdateDropLocationEvent, location);
+
+            // Reply that we are ready to copy the dragged data
+            XEvent reply;
+            memset(&reply, 0, sizeof(reply));
+
+            reply.type = ClientMessage;
+            reply.xclient.window       = event->xclient.data.l[0];
+            reply.xclient.message_type = this->XdndStatusAtom;
+            reply.xclient.format       = 32;
+            reply.xclient.data.l[0]    = this->WindowId;
+            reply.xclient.data.l[1]    = 1; // Always accept the dnd with no rectangle
+            reply.xclient.data.l[2]    = 0; // Specify an empty rectangle
+            reply.xclient.data.l[3]    = 0;
+            reply.xclient.data.l[4]    = this->XdndActionCopyAtom;
+
+            XSendEvent(this->DisplayId, event->xclient.data.l[0], False, NoEventMask, &reply);
+            XFlush(this->DisplayId);
+        }
+        else if (event->xclient.message_type == this->XdndDropAtom)
+        {
+            // Item dropped in the window
+            // Store the source of the drag and drop
+            this->XdndSource = event->xclient.data.l[0];
 
-        reply.type = ClientMessage;
-        reply.xclient.window = event->xclient.data.l[0];
-        reply.xclient.message_type = this->XdndStatusAtom;
-        reply.xclient.format = 32;
-        reply.xclient.data.l[0] = this->WindowId;
-        reply.xclient.data.l[1] = 1; // Always accept the dnd with no rectangle
-        reply.xclient.data.l[2] = 0; // Specify an empty rectangle
-        reply.xclient.data.l[3] = 0;
-        reply.xclient.data.l[4] = this->XdndActionCopyAtom;
-
-        XSendEvent(this->DisplayId, event->xclient.data.l[0], False, NoEventMask, &reply);
-        XFlush(this->DisplayId);
-      }
-      else if (event->xclient.message_type == this->XdndDropAtom)
-      {
-        // Item dropped in the window
-        // Store the source of the drag and drop
-        this->XdndSource = event->xclient.data.l[0];
-
-        // Ask for a conversion of the selection. This will trigger a SelectioNotify event later.
-        Atom xdndSelectionAtom = XInternAtom(this->DisplayId, "XdndSelection", False);
-        XConvertSelection(this->DisplayId, xdndSelectionAtom,
+            // Ask for a conversion of the selection. This will trigger a SelectioNotify event later.
+            Atom xdndSelectionAtom = XInternAtom(this->DisplayId, "XdndSelection", False);
+            XConvertSelection(this->DisplayId, xdndSelectionAtom,
           XInternAtom(this->DisplayId, "UTF8_STRING", False), xdndSelectionAtom, this->WindowId,
           CurrentTime);
-      }
-      else if (static_cast<Atom>(event->xclient.data.l[0]) == this->KillAtom)
-      {
-        this->ExitCallback();
-      }
+        }
+        else if (static_cast<Atom>(event->xclient.data.l[0]) == this->KillAtom)
+        {
+            this->ExitCallback();
+        }
     }
     break;
-  }
+    }
 }
 
 //------------------------------------------------------------------------------
-void vtkXRenderWindowInteractor2::GetMousePosition(int* x, int* y)
+void
+vtkXRenderWindowInteractor2::GetMousePosition(int* x, int* y)
 {
-  Window root, child;
-  int root_x, root_y;
-  unsigned int keys;
+    Window       root, child;
+    int          root_x, root_y;
+    unsigned int keys;
 
-  XQueryPointer(this->DisplayId, this->WindowId, &root, &child, &root_x, &root_y, x, y, &keys);
+    XQueryPointer(this->DisplayId, this->WindowId, &root, &child, &root_x, &root_y, x, y, &keys);
 
-  *y = this->Size[1] - *y - 1;
+    *y = this->Size[1] - *y - 1;
 }
 
 //------------------------------------------------------------------------------
diff --git a/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.h b/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.h
index 3afa9f759518dd375b498a828050f6a45b9304a3..026ea6ca1800c40bf10c02da1eadb95ed441d0b2 100644
--- a/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.h
+++ b/Source/SimulationManager/VTKRenderer/imstkVtkXRenderWindowInteractor2.h
@@ -40,102 +40,102 @@ class vtkXRenderWindowInteractor2Internals;
 class vtkXRenderWindowInteractor2 : public vtkRenderWindowInteractor
 {
 public:
-  static vtkXRenderWindowInteractor2* New();
-  vtkTypeMacro(vtkXRenderWindowInteractor2, vtkRenderWindowInteractor);
-  void PrintSelf(ostream& os, vtkIndent indent) override;
-
-  /**
-   * Initializes the event handlers without an XtAppContext.  This is
-   * good for when you don't have a user interface, but you still
-   * want to have mouse interaction.
-   */
-  void Initialize() override;
-
-  /**
-   * Break the event loop on 'q','e' keypress. Want more ???
-   */
-  void TerminateApp() override;
-
-  /**
-   * Run the event loop and return. This is provided so that you can
-   * implement your own event loop but yet use the vtk event handling as
-   * well.
-   */
-  void ProcessEvents() override;
-
-  //@{
-  /**
-   * Enable/Disable interactions.  By default interactors are enabled when
-   * initialized.  Initialize() must be called prior to enabling/disabling
-   * interaction. These methods are used when a window/widget is being
-   * shared by multiple renderers and interactors.  This allows a "modal"
-   * display where one interactor is active when its data is to be displayed
-   * and all other interactors associated with the widget are disabled
-   * when their data is not displayed.
-   */
-  void Enable() override;
-  void Disable() override;
-  //@}
-
-  /**
-   * Update the Size data member and set the associated RenderWindow's
-   * size.
-   */
-  void UpdateSize(int, int) override;
-
-  /**
-   * Re-defines virtual function to get mouse position by querying X-server.
-   */
-  void GetMousePosition(int* x, int* y) override;
-
-  void DispatchEvent(XEvent*);
+    static vtkXRenderWindowInteractor2* New();
+    vtkTypeMacro(vtkXRenderWindowInteractor2, vtkRenderWindowInteractor);
+    void PrintSelf(ostream& os, vtkIndent indent) override;
+
+    /**
+     * Initializes the event handlers without an XtAppContext.  This is
+     * good for when you don't have a user interface, but you still
+     * want to have mouse interaction.
+     */
+    void Initialize() override;
+
+    /**
+     * Break the event loop on 'q','e' keypress. Want more ???
+     */
+    void TerminateApp() override;
+
+    /**
+     * Run the event loop and return. This is provided so that you can
+     * implement your own event loop but yet use the vtk event handling as
+     * well.
+     */
+    void ProcessEvents() override;
+
+    //@{
+    /**
+     * Enable/Disable interactions.  By default interactors are enabled when
+     * initialized.  Initialize() must be called prior to enabling/disabling
+     * interaction. These methods are used when a window/widget is being
+     * shared by multiple renderers and interactors.  This allows a "modal"
+     * display where one interactor is active when its data is to be displayed
+     * and all other interactors associated with the widget are disabled
+     * when their data is not displayed.
+     */
+    void Enable() override;
+    void Disable() override;
+    //@}
+
+    /**
+     * Update the Size data member and set the associated RenderWindow's
+     * size.
+     */
+    void UpdateSize(int, int) override;
+
+    /**
+     * Re-defines virtual function to get mouse position by querying X-server.
+     */
+    void GetMousePosition(int* x, int* y) override;
+
+    void DispatchEvent(XEvent*);
 
 protected:
-  vtkXRenderWindowInteractor2();
-  ~vtkXRenderWindowInteractor2() override;
-
-  /**
-   * Update the Size data member and set the associated RenderWindow's
-   * size but do not resize the XWindow.
-   */
-  void UpdateSizeNoXResize(int, int);
-
-  // Using static here to avoid destroying context when many apps are open:
-  static int NumAppInitialized;
-
-  Display* DisplayId;
-  Window WindowId;
-  Atom KillAtom;
-  int PositionBeforeStereo[2];
-  vtkXRenderWindowInteractor2Internals* Internal;
-
-  // Drag and drop related
-  Window XdndSource;
-  Atom XdndPositionAtom;
-  Atom XdndDropAtom;
-  Atom XdndActionCopyAtom;
-  Atom XdndStatusAtom;
-  Atom XdndFinishedAtom;
-
-  //@{
-  /**
-   * X-specific internal timer methods. See the superclass for detailed
-   * documentation.
-   */
-  int InternalCreateTimer(int timerId, int timerType, unsigned long duration) override;
-  int InternalDestroyTimer(int platformTimerId) override;
-  //@}
-
-  void FireTimers();
-
-  /**
-   * This will start up the X event loop and never return. If you
-   * call this method it will loop processing X events until the
-   * application is exited.
-   */
-  void StartEventLoop() override;
+    vtkXRenderWindowInteractor2();
+    ~vtkXRenderWindowInteractor2() override;
+
+    /**
+     * Update the Size data member and set the associated RenderWindow's
+     * size but do not resize the XWindow.
+     */
+    void UpdateSizeNoXResize(int, int);
+
+    // Using static here to avoid destroying context when many apps are open:
+    static int NumAppInitialized;
+
+    Display* DisplayId;
+    Window   WindowId;
+    Atom     KillAtom;
+    int      PositionBeforeStereo[2];
+    vtkXRenderWindowInteractor2Internals* Internal;
+
+    // Drag and drop related
+    Window XdndSource;
+    Atom   XdndPositionAtom;
+    Atom   XdndDropAtom;
+    Atom   XdndActionCopyAtom;
+    Atom   XdndStatusAtom;
+    Atom   XdndFinishedAtom;
+
+    //@{
+    /**
+     * X-specific internal timer methods. See the superclass for detailed
+     * documentation.
+     */
+    int InternalCreateTimer(int timerId, int timerType, unsigned long duration) override;
+    int InternalDestroyTimer(int platformTimerId) override;
+    //@}
+
+    void FireTimers();
+
+    /**
+     * This will start up the X event loop and never return. If you
+     * call this method it will loop processing X events until the
+     * application is exited.
+     */
+    void StartEventLoop() override;
 
 private:
-  vtkXRenderWindowInteractor2(const vtkXRenderWindowInteractor2&) = delete;
-  void operator=(const vtkXRenderWindowInteractor2&) = delete;
+    vtkXRenderWindowInteractor2(const vtkXRenderWindowInteractor2&) = delete;
+    void operator=(const vtkXRenderWindowInteractor2&) = delete;
 };
\ No newline at end of file
diff --git a/Source/SimulationManager/imstkSimulationManager.h b/Source/SimulationManager/imstkSimulationManager.h
index b2949b425e64667821bca2fbd63a654e36a40b69..b8a1a364a545e6e2bf3a3bc892a1c05bab7e1f1a 100644
--- a/Source/SimulationManager/imstkSimulationManager.h
+++ b/Source/SimulationManager/imstkSimulationManager.h
@@ -114,12 +114,12 @@ protected:
 protected:
     std::vector<std::shared_ptr<Viewer>> m_viewers;
 
+    std::unordered_map<Module*,bool> m_running;
+
     std::vector<std::shared_ptr<Module>> m_syncModules;      ///> Modules called once per update
     std::vector<std::shared_ptr<Module>> m_asyncModules;     ///> Modules that run on completely other threads without restraint
     std::vector<std::shared_ptr<Module>> m_adaptiveModules;  ///> Modules that update adpatively to keep up with real time
 
-    std::unordered_map<Module*,bool> m_running;
-
     ThreadingType m_threadType = ThreadingType::STL;
     double m_desiredDt = 0.003; // Desired timestep
     double m_dt       = 0.0;    // Actual timestep