diff --git a/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp b/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp
index 8aeeca5ed7163839a2887c5e84034b93416a3f79..7c1253ea4e349f5c6783632a7d4629acfc0fe123 100644
--- a/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp
+++ b/Examples/ComputeGraph/Benchmarking/computeGraphBenchmarkExample.cpp
@@ -98,9 +98,9 @@ main()
         sims[i].params->enableConstraint(PbdConstraint::Type::Bend, static_cast<double>(i) * 0.1 / numStrings + 0.001);
         sims[i].params->m_fixedNodeIds     = { 0 };
         sims[i].params->m_uniformMassValue = 5.0;
-        sims[i].params->m_gravity   = Vec3d(0, -9.8, 0);
-        sims[i].params->m_defaultDt = 0.0005;
-        sims[i].params->m_iterations   = 5;
+        sims[i].params->m_gravity    = Vec3d(0, -9.8, 0);
+        sims[i].params->m_defaultDt  = 0.0005;
+        sims[i].params->m_iterations = 5;
 
         // Set the parameters
         sims[i].model->configure(sims[i].params);
@@ -129,18 +129,18 @@ main()
     const double radius     = 1.5;
     auto         movePoints =
         [&sims, &t, dt, radius](Module* module)
-    {
-        for (unsigned int i = 0; i < sims.size(); i++)
         {
-            Vec3d pos = sims[i].model->getCurrentState()->getVertexPosition(0);
-            // Move in circle, derivatives of parametric eq of circle
-            sims[i].model->getCurrentState()->setVertexPosition(0, imstk::Vec3d(
+            for (unsigned int i = 0; i < sims.size(); i++)
+            {
+                Vec3d pos = sims[i].model->getCurrentState()->getVertexPosition(0);
+                // Move in circle, derivatives of parametric eq of circle
+                sims[i].model->getCurrentState()->setVertexPosition(0, imstk::Vec3d(
                 pos.x() + -std::sin(t) * radius * dt,
                 pos.y(),
                 pos.z() + std::cos(t) * radius * dt));
-        }
-        t += dt;
-    };
+            }
+            t += dt;
+        };
     simManager->getSceneManager(scene)->setPostUpdateCallback(movePoints);
 
     // Start
diff --git a/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp b/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp
index 65e03f79d75b152d024e615ba5556b7b22611aad..a4737aa6ae15fefba536540c46e811c81d043782 100644
--- a/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp
+++ b/Examples/ComputeGraph/Configuration/computeGraphConfigureExample.cpp
@@ -38,13 +38,14 @@
 
 using namespace imstk;
 
-static std::unique_ptr<SurfaceMesh> makeCloth(
+static std::unique_ptr<SurfaceMesh>
+makeCloth(
     const double width, const double height, const int nRows, const int nCols)
 {
     // Create surface mesh
     std::unique_ptr<SurfaceMesh> clothMesh = std::make_unique<SurfaceMesh>();
-    StdVectorOfVec3d vertList;
-    
+    StdVectorOfVec3d             vertList;
+
     vertList.resize(nRows * nCols);
     const double dy = width / (double)(nCols - 1);
     const double dx = height / (double)(nRows - 1);
@@ -91,7 +92,8 @@ static std::unique_ptr<SurfaceMesh> makeCloth(
     return clothMesh;
 }
 
-static std::shared_ptr<PbdObject> makeClothObj(const std::string& name, double width, double height, int nRows, int nCols)
+static std::shared_ptr<PbdObject>
+makeClothObj(const std::string& name, double width, double height, int nRows, int nCols)
 {
     auto clothObj = std::make_shared<PbdObject>(name);
 
@@ -100,9 +102,9 @@ static std::shared_ptr<PbdObject> makeClothObj(const std::string& name, double w
     // Setup the Parameters
     auto pbdParams = std::make_shared<PBDModelConfig>();
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -9.8, 0);
-    pbdParams->m_defaultDt = 0.005;
-    pbdParams->m_iterations   = 5;
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
+    pbdParams->m_defaultDt  = 0.005;
+    pbdParams->m_iterations = 5;
     pbdParams->m_viscousDampingCoeff = 0.05;
     pbdParams->enableConstraint(PbdConstraint::Type::Distance, 0.1);
     pbdParams->enableConstraint(PbdConstraint::Type::Dihedral, 0.001);
@@ -144,10 +146,10 @@ main()
     auto simManager = std::make_shared<SimulationManager>();
     auto scene      = simManager->createNewScene("PBDCloth");
 
-    const double width  = 10.0;
-    const double height = 10.0;
-    const int nRows = 16;
-    const int nCols = 16;
+    const double               width    = 10.0;
+    const double               height   = 10.0;
+    const int                  nRows    = 16;
+    const int                  nCols    = 16;
     std::shared_ptr<PbdObject> clothObj = makeClothObj("Cloth", width, height, nRows, nCols);
     scene->addSceneObject(clothObj);
 
@@ -174,33 +176,32 @@ main()
     {
         std::shared_ptr<PbdModel> pbdModel = clothObj->getPbdModel();
         scene->setComputeGraphConfigureCallback([&](Scene* scene)
-            {
-                // Get the graph
-                std::shared_ptr<ComputeGraph> graph = scene->getComputeGraph();
-
-                // First write the graph before we make modifications, just to show the changes
-                ComputeGraphVizWriter writer;
-                writer.setInput(graph);
-                writer.setFileName("computeGraphConfigureExampleOld.svg");
-                writer.write();
-
-                std::shared_ptr<ComputeNode> printVelocities = std::make_shared<ComputeNode>([&]()
-                    {
-                        // Make the timestep a function of max velocity
-                        const StdVectorOfVec3d& velocities = *pbdModel->getCurrentState()->getVelocities();
-                        for (size_t i = 0; i < velocities.size(); i++)
-                        {
-                            printf("Velocity: %f, %f, %f\n", velocities[i].x(), velocities[i].y(), velocities[i].z());
-                        }
+        {
+            // Get the graph
+            std::shared_ptr<ComputeGraph> graph = scene->getComputeGraph();
 
+            // First write the graph before we make modifications, just to show the changes
+            ComputeGraphVizWriter writer;
+            writer.setInput(graph);
+            writer.setFileName("computeGraphConfigureExampleOld.svg");
+            writer.write();
+
+            std::shared_ptr<ComputeNode> printVelocities = std::make_shared<ComputeNode>([&]()
+            {
+                // Make the timestep a function of max velocity
+                const StdVectorOfVec3d& velocities = *pbdModel->getCurrentState()->getVelocities();
+                for (size_t i = 0; i < velocities.size(); i++)
+                {
+                    printf("Velocity: %f, %f, %f\n", velocities[i].x(), velocities[i].y(), velocities[i].z());
+                }
                     }, "PrintVelocities");
 
-                // After IntegratePosition
-                graph->insertAfter(pbdModel->getIntegratePositionNode(), printVelocities);
+            // After IntegratePosition
+            graph->insertAfter(pbdModel->getIntegratePositionNode(), printVelocities);
 
-                // Write the modified graph
-                writer.setFileName("computeGraphConfigureExampleNew.svg");
-                writer.write();
+            // Write the modified graph
+            writer.setFileName("computeGraphConfigureExampleNew.svg");
+            writer.write();
             });
     }
 
diff --git a/Examples/DeformableBody/DeformableBodyExample.cpp b/Examples/DeformableBody/DeformableBodyExample.cpp
index 1ba324f17a5d4ea41068b21988ec33eb888a42ef..86e6e0ad4b082736a77924f154c87e356eb41b9b 100644
--- a/Examples/DeformableBody/DeformableBodyExample.cpp
+++ b/Examples/DeformableBody/DeformableBodyExample.cpp
@@ -105,30 +105,6 @@ main()
     planeObj->setCollidingGeometry(planeGeom);
     scene->addSceneObject(planeObj);
 
-    // create a nonlinear system
-    auto nlSystem = std::make_shared<NonLinearSystem<SparseMatrixd>>(
-        dynaModel->getFunction(),
-        dynaModel->getFunctionGradient());
-
-    nlSystem->setUnknownVector(dynaModel->getUnknownVec());
-    nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
-    nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction());
-
-    // create a linear solver
-    auto linSolver = std::make_shared<ConjugateGradient>();
-
-    if (linSolver->getType() == imstk::LinearSolver<imstk::SparseMatrixd>::Type::GaussSeidel
-        && dynaModel->isFixedBCImplemented())
-    {
-        LOG(WARNING) << "The GS solver may not be viable!";
-    }
-
-    // create a non-linear solver and add to the scene
-    auto nlSolver = std::make_shared<NewtonSolver<SparseMatrixd>>();
-    nlSolver->setLinearSolver(linSolver);
-    nlSolver->setSystem(nlSystem);
-    scene->addNonlinearSolver(nlSolver);
-
     // Light
     auto light = std::make_shared<DirectionalLight>("light");
     light->setFocalPoint(Vec3d(5, -8, -5));
diff --git a/Examples/MultipleScenes/multipleScenes.cpp b/Examples/MultipleScenes/multipleScenes.cpp
index fbe49abf21bf16c1b481953b7a655d80f876ec7a..e120419ee22a37d9f17b185365812f075b0f0217 100644
--- a/Examples/MultipleScenes/multipleScenes.cpp
+++ b/Examples/MultipleScenes/multipleScenes.cpp
@@ -72,7 +72,7 @@ createSoftBodyScene(std::shared_ptr<SimulationManager> simManager, const char* s
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity = Vec3d(0, -9.8, 0);
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
     pbdParams->m_iterations = 45;
 
     // Set the parameters
@@ -163,9 +163,9 @@ createClothScene(std::shared_ptr<SimulationManager> simManager, const char* scen
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -9.8, 0);
-    pbdParams->m_defaultDt = 0.03;
-    pbdParams->m_iterations   = 5;
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
+    pbdParams->m_defaultDt  = 0.03;
+    pbdParams->m_iterations = 5;
 
     // Set the parameters
     pbdModel->configure(pbdParams);
diff --git a/Examples/NonlinearSolver/nonlinearSolver.cpp b/Examples/NonlinearSolver/nonlinearSolver.cpp
index 1725f74c8351bfb07f27fed3fbb7c273434782d9..1e50fa5a3f7ba643b2ef9eb0466e52ce2ab20dd6 100644
--- a/Examples/NonlinearSolver/nonlinearSolver.cpp
+++ b/Examples/NonlinearSolver/nonlinearSolver.cpp
@@ -34,40 +34,40 @@ using namespace imstk;
 int
 main(int argc, char** argv)
 {
-    const int N = 2;
-    auto      x = Vectord(N);
+    const int N  = 2;
+    auto      x  = Vectord(N);
     auto      xe = Vectord(N);
-    auto      y = Vectord(N);
-    auto      A = Matrixd(N, N);
+    auto      y  = Vectord(N);
+    auto      A  = Matrixd(N, N);
 
-    x[0] = 100.0;
-    x[1] = 100.0;
+    x[0]  = 100.0;
+    x[1]  = 100.0;
     xe[0] = 1.0;
     xe[1] = 10.0;
 
     auto func = [&y](const Vectord& x, const bool isSemiImplicit) -> const Vectord& {
-        // auto y = Vectord(x.size());
-        y[0] = x[0] * x[0] - 1.0;
-        y[1] = x[1] * x[1] - 100.0;
+                    // auto y = Vectord(x.size());
+                    y[0] = x[0] * x[0] - 1.0;
+                    y[1] = x[1] * x[1] - 100.0;
 
-        return y;
-    };
+                    return y;
+                };
 
     auto jac = [&A](const Vectord& x) -> const Matrixd& {
-        // auto A = Matrixd(x.size(), x.size());
-        A(0, 0) = 2 * x[0];
-        A(0, 1) = 0.0;
-        A(1, 0) = 0.0;
-        A(1, 1) = 2 * x[1];
+                   // auto A = Matrixd(x.size(), x.size());
+                   A(0, 0) = 2 * x[0];
+                   A(0, 1) = 0.0;
+                   A(1, 0) = 0.0;
+                   A(1, 1) = 2 * x[1];
 
-        return A;
-    };
+                   return A;
+               };
 
     auto updateX = [&x](const Vectord& du, const bool isSemiImplicit)
-    {
-        x -= du;
-        return;
-    };
+                   {
+                       x -= du;
+                       return;
+                   };
 
     auto updateXold = [](void) {};
 
@@ -84,8 +84,8 @@ main(int argc, char** argv)
     nlSolver->setSystem(nlSystem);
     nlSolver->setLinearSolver(linSolver);
 
-    std::cout << "init_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl;
+    std::cout << "init_error = " << std::setprecision(12) << std::scientific << (x - xe).norm() << std::endl;
     nlSolver->solve();
 
-    std::cout << "final_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl;
+    std::cout << "final_error = " << std::setprecision(12) << std::scientific << (x - xe).norm() << std::endl;
 }
diff --git a/Examples/PBD/PBDCloth/pbdClothExample.cpp b/Examples/PBD/PBDCloth/pbdClothExample.cpp
index 0477de14088ebff2dd143e1e77ff407c877f9182..c76b56a558a1d6fb254781d2d0962a357bd3663c 100644
--- a/Examples/PBD/PBDCloth/pbdClothExample.cpp
+++ b/Examples/PBD/PBDCloth/pbdClothExample.cpp
@@ -110,9 +110,9 @@ main()
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -9.8, 0);
-    pbdParams->m_defaultDt = 0.005;
-    pbdParams->m_iterations   = 5;
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
+    pbdParams->m_defaultDt  = 0.005;
+    pbdParams->m_iterations = 5;
 
     // Set the parameters
     pbdModel->configure(pbdParams);
diff --git a/Examples/PBD/PBDCollisionMultipleObjects/PBDCollisionMultipleObjectsExample.cpp b/Examples/PBD/PBDCollisionMultipleObjects/PBDCollisionMultipleObjectsExample.cpp
index 98da44f06aa296e1034e8f41faf58d981de1e9f3..c6925ef6d4c53748969f08d6115e035f807e4a3d 100644
--- a/Examples/PBD/PBDCollisionMultipleObjects/PBDCollisionMultipleObjectsExample.cpp
+++ b/Examples/PBD/PBDCollisionMultipleObjects/PBDCollisionMultipleObjectsExample.cpp
@@ -125,8 +125,8 @@ generateDragon(const std::shared_ptr<imstk::Scene>& scene,
 
     // Other parameters
     pbdParams->m_uniformMassValue = 5.0;
-    pbdParams->m_gravity   = Vec3d(0, -1.0, 0);
-    pbdParams->m_defaultDt = 0.01;
+    pbdParams->m_gravity    = Vec3d(0, -1.0, 0);
+    pbdParams->m_defaultDt  = 0.01;
     pbdParams->m_iterations = 20;
     pbdParams->collisionParams->m_proximity = 0.5;
 
@@ -222,7 +222,7 @@ main()
     const double distanceY  = 5.0;
     const double minHeight  = -5.0;
 
-    std::vector<std::shared_ptr<PbdObject>>   pbdObjs;
+    std::vector<std::shared_ptr<PbdObject>> pbdObjs;
 
 #ifdef BIG_SCENE
     for (int i = -expandsXZ; i < expandsXZ; ++i)
diff --git a/Examples/PBD/PBDCollisionOneObject/PBDCollisionOneObjectExample.cpp b/Examples/PBD/PBDCollisionOneObject/PBDCollisionOneObjectExample.cpp
index 8616bf08bbbb6099dc4b34f6d6b9afeb67e1f6df..6fc7c41ba69ebe94343a148c7e0d89bec49dee80 100644
--- a/Examples/PBD/PBDCollisionOneObject/PBDCollisionOneObjectExample.cpp
+++ b/Examples/PBD/PBDCollisionOneObject/PBDCollisionOneObjectExample.cpp
@@ -85,9 +85,9 @@ main()
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -10.0, 0);
-    pbdParams->m_defaultDt = 0.01;
-    pbdParams->m_iterations   = 5;
+    pbdParams->m_gravity    = Vec3d(0, -10.0, 0);
+    pbdParams->m_defaultDt  = 0.01;
+    pbdParams->m_iterations = 5;
     pbdParams->collisionParams->m_proximity = 0.3;
     pbdParams->collisionParams->m_stiffness = 0.1;
 
@@ -149,7 +149,7 @@ main()
     // configure model
     auto pbdParams2 = std::make_shared<PBDModelConfig>();
     pbdParams2->m_uniformMassValue = 0.0;
-    pbdParams2->m_iterations = 0;
+    pbdParams2->m_iterations       = 0;
     pbdParams2->collisionParams->m_proximity = -0.1;
 
     // Set the parameters
diff --git a/Examples/PBD/PBDCollisionStairs/PBDCollisionStairsExample.cpp b/Examples/PBD/PBDCollisionStairs/PBDCollisionStairsExample.cpp
index bbc2c685a17dfe72de8b4d455aad6c9b1b8afaa6..ef8da4689b69190b0f8f44169227602646ae0c72 100644
--- a/Examples/PBD/PBDCollisionStairs/PBDCollisionStairsExample.cpp
+++ b/Examples/PBD/PBDCollisionStairs/PBDCollisionStairsExample.cpp
@@ -113,8 +113,8 @@ makeDragonPbdObject(const std::string& name)
     pbdParams->enableFEMConstraint(PbdConstraint::Type::FEMTet,
         PbdFEMConstraint::MaterialType::StVK);
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -10.0, 0);
-    pbdParams->m_defaultDt = 0.01;
+    pbdParams->m_gravity    = Vec3d(0, -10.0, 0);
+    pbdParams->m_defaultDt  = 0.01;
     pbdParams->m_iterations = 10;
     pbdParams->collisionParams->m_proximity = 0.3;
     pbdParams->collisionParams->m_stiffness = 0.1;
diff --git a/Examples/PBD/PBDDeformableObject/PBD3DDeformableObject.cpp b/Examples/PBD/PBDDeformableObject/PBD3DDeformableObject.cpp
index b52a57a9e1a568589adf1847e31d0e61cfb20180..b07e73792aabf9ecef7a0c46f7511e5f9b02ceec 100644
--- a/Examples/PBD/PBDDeformableObject/PBD3DDeformableObject.cpp
+++ b/Examples/PBD/PBDDeformableObject/PBD3DDeformableObject.cpp
@@ -73,7 +73,7 @@ main()
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity = Vec3d(0, -9.8, 0);
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
     pbdParams->m_iterations = 45;
 
     // Set the parameters
diff --git a/Examples/PBD/PBDFluids/PBDFluidsExample.cpp b/Examples/PBD/PBDFluids/PBDFluidsExample.cpp
index 46c002e16c40e76d6be38ab4cab905fe98bb0195..002a5b4e255e426ecfaf5f1908609613b2680178 100644
--- a/Examples/PBD/PBDFluids/PBDFluidsExample.cpp
+++ b/Examples/PBD/PBDFluids/PBDFluidsExample.cpp
@@ -76,9 +76,9 @@ main()
 
     // Other parameters
     pbdParams->m_uniformMassValue = 1.0;
-    pbdParams->m_gravity   = Vec3d(0, -9.8, 0);
-    pbdParams->m_defaultDt = 0.01;
-    pbdParams->m_iterations   = 2;
+    pbdParams->m_gravity    = Vec3d(0, -9.8, 0);
+    pbdParams->m_defaultDt  = 0.01;
+    pbdParams->m_iterations = 2;
     pbdParams->collisionParams->m_proximity = 0.01;
 
     // Set the parameters
diff --git a/Examples/PBD/PBDString/pbdStringExample.cpp b/Examples/PBD/PBDString/pbdStringExample.cpp
index 9d34aa3120c4d61aea3d14932091052d1ece55cf..b3b17a9eafb256b3461d54bf6420da5a81d54692 100644
--- a/Examples/PBD/PBDString/pbdStringExample.cpp
+++ b/Examples/PBD/PBDString/pbdStringExample.cpp
@@ -95,9 +95,9 @@ main()
         sims[i].params->enableConstraint(PbdConstraint::Type::Bend, static_cast<double>(i) * 0.1 / numStrings + 0.001);
         sims[i].params->m_fixedNodeIds     = { 0 };
         sims[i].params->m_uniformMassValue = 5.0;
-        sims[i].params->m_gravity   = Vec3d(0, -9.8, 0);
-        sims[i].params->m_defaultDt = 0.0005;
-        sims[i].params->m_iterations   = 5;
+        sims[i].params->m_gravity    = Vec3d(0, -9.8, 0);
+        sims[i].params->m_defaultDt  = 0.0005;
+        sims[i].params->m_iterations = 5;
 
         // Set the parameters
         sims[i].model->configure(sims[i].params);
diff --git a/Source/CollisionDetection/CollisionDetection/imstkCollisionDetection.h b/Source/CollisionDetection/CollisionDetection/imstkCollisionDetection.h
index b7481a5ae7712c647f5d5ce2b3c9c6280b049ab4..bbaebce35387da5092dfa3f5682a870862d63b49 100644
--- a/Source/CollisionDetection/CollisionDetection/imstkCollisionDetection.h
+++ b/Source/CollisionDetection/CollisionDetection/imstkCollisionDetection.h
@@ -64,7 +64,7 @@ public:
         BidirectionalPlaneToSphere,
         SphereToCylinder,
         SphereToSphere,
-        
+
         // Image based CD
         SignedDistanceField,
 
diff --git a/Source/CollisionDetection/imstkCDObjectFactory.cpp b/Source/CollisionDetection/imstkCDObjectFactory.cpp
index b31e5c5987719710d0787a0f4738bcead82b4882..668990d220005b1e0974f328ba32820441ed2c1e 100644
--- a/Source/CollisionDetection/imstkCDObjectFactory.cpp
+++ b/Source/CollisionDetection/imstkCDObjectFactory.cpp
@@ -61,33 +61,33 @@ namespace imstk
 {
 std::shared_ptr<CollisionDetection>
 makeCollisionDetectionObject(const CollisionDetection::Type type,
-    std::shared_ptr<Geometry> collidingGeometryA,
-    std::shared_ptr<Geometry> collidingGeometryB,
-    std::shared_ptr<CollisionData> colData)
+                             std::shared_ptr<Geometry>      collidingGeometryA,
+                             std::shared_ptr<Geometry>      collidingGeometryB,
+                             std::shared_ptr<CollisionData> colData)
 {
     switch (type)
     {
-        // Points to objects
+    // Points to objects
     case CollisionDetection::Type::PointSetToSphere:
     {
         auto pointset = std::dynamic_pointer_cast<PointSet>(collidingGeometryA);
         auto sphere   = std::dynamic_pointer_cast<Sphere>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(pointset, sphere)
-            return std::make_shared<PointSetToSphereCD>(pointset, sphere, colData);
+        return std::make_shared<PointSetToSphereCD>(pointset, sphere, colData);
     }
     case CollisionDetection::Type::PointSetToPlane:
     {
         auto pointset = std::dynamic_pointer_cast<PointSet>(collidingGeometryA);
         auto plane    = std::dynamic_pointer_cast<Plane>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(pointset, plane)
-            return std::make_shared<PointSetToPlaneCD>(pointset, plane, colData);
+        return std::make_shared<PointSetToPlaneCD>(pointset, plane, colData);
     }
     case CollisionDetection::Type::PointSetToCapsule:
     {
         auto pointset = std::dynamic_pointer_cast<PointSet>(collidingGeometryA);
         auto capsule  = std::dynamic_pointer_cast<Capsule>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(pointset, capsule)
-            return std::make_shared<PointSetToCapsuleCD>(pointset, capsule, colData);
+        return std::make_shared<PointSetToCapsuleCD>(pointset, capsule, colData);
     }
 #if 0
     case CollisionDetection::ype::PointSetToSpherePicking:
@@ -95,24 +95,24 @@ makeCollisionDetectionObject(const CollisionDetection::Type type,
         auto pointset = std::dynamic_pointer_cast<PointSet>(objA->getCollidingGeometry());
         auto sphere   = std::dynamic_pointer_cast<Sphere>(objB->getCollidingGeometry());
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(pointset, sphere)
-            return std::make_shared<PointSetToSpherePickingCD>(pointset, sphere, colData);
+        return std::make_shared<PointSetToSpherePickingCD>(pointset, sphere, colData);
     }
 #endif
     case CollisionDetection::Type::PointSetToSurfaceMesh:
     {
-        auto        pointset = std::dynamic_pointer_cast<PointSet>(collidingGeometryA);
-        auto        triMesh  = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
+        auto pointset = std::dynamic_pointer_cast<PointSet>(collidingGeometryA);
+        auto triMesh  = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(pointset, triMesh)
-            CollisionDetection::addCollisionPairToOctree(collidingGeometryA, collidingGeometryB, type, colData);
+        CollisionDetection::addCollisionPairToOctree(collidingGeometryA, collidingGeometryB, type, colData);
         return std::make_shared<PointSetToSurfaceMeshCD>(pointset, triMesh, colData);
     }
     // Mesh to mesh
     case CollisionDetection::Type::SurfaceMeshToSurfaceMesh:
     {
-        auto        meshA = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryA);
-        auto        meshB = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
+        auto meshA = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryA);
+        auto meshB = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(meshA, meshB)
-            CollisionDetection::addCollisionPairToOctree(collidingGeometryA, collidingGeometryB, type, colData);
+        CollisionDetection::addCollisionPairToOctree(collidingGeometryA, collidingGeometryB, type, colData);
         return std::make_shared<SurfaceMeshToSurfaceMeshCD>(meshA, meshB, colData);
     }
     case CollisionDetection::Type::SurfaceMeshToSurfaceMeshCCD:
@@ -120,7 +120,7 @@ makeCollisionDetectionObject(const CollisionDetection::Type type,
         auto meshA = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryA);
         auto meshB = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(meshA, meshB)
-            return std::make_shared<SurfaceMeshToSurfaceMeshCCD>(meshA, meshB, colData);
+        return std::make_shared<SurfaceMeshToSurfaceMeshCCD>(meshA, meshB, colData);
     }
 
     case CollisionDetection::Type::VolumeMeshToVolumeMesh:
@@ -128,7 +128,7 @@ makeCollisionDetectionObject(const CollisionDetection::Type type,
         auto tet1 = std::dynamic_pointer_cast<TetrahedralMesh>(collidingGeometryA);
         auto tet2 = std::dynamic_pointer_cast<TetrahedralMesh>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(tet1, tet2)
-            return std::make_shared<TetraToTetraCD>(tet1, tet2, colData);
+        return std::make_shared<TetraToTetraCD>(tet1, tet2, colData);
     }
 
     // Analytical object to analytical object
@@ -137,34 +137,34 @@ makeCollisionDetectionObject(const CollisionDetection::Type type,
         auto plane  = std::dynamic_pointer_cast<Plane>(collidingGeometryA);
         auto sphere = std::dynamic_pointer_cast<Sphere>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(plane, sphere)
-            return std::make_shared<UnidirectionalPlaneToSphereCD>(plane, sphere, colData);
+        return std::make_shared<UnidirectionalPlaneToSphereCD>(plane, sphere, colData);
     }
     case CollisionDetection::Type::BidirectionalPlaneToSphere:
     {
         auto plane  = std::dynamic_pointer_cast<Plane>(collidingGeometryA);
         auto sphere = std::dynamic_pointer_cast<Sphere>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(plane, sphere)
-            return std::make_shared<BidirectionalPlaneToSphereCD>(plane, sphere, colData);
+        return std::make_shared<BidirectionalPlaneToSphereCD>(plane, sphere, colData);
     }
     case CollisionDetection::Type::SphereToSphere:
     {
         auto sphereA = std::dynamic_pointer_cast<Sphere>(collidingGeometryA);
         auto sphereB = std::dynamic_pointer_cast<Sphere>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(sphereA, sphereB)
-            return std::make_shared<SphereToSphereCD>(sphereA, sphereB, colData);
+        return std::make_shared<SphereToSphereCD>(sphereA, sphereB, colData);
     }
     case CollisionDetection::Type::SphereToCylinder:
     {
         auto sphere   = std::dynamic_pointer_cast<Sphere>(collidingGeometryB);
         auto cylinder = std::dynamic_pointer_cast<Cylinder>(collidingGeometryA);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(sphere, cylinder)
-            return std::make_shared<SphereToCylinderCD>(sphere, cylinder, colData);
+        return std::make_shared<SphereToCylinderCD>(sphere, cylinder, colData);
     }
     case CollisionDetection::Type::MeshToMeshBruteForce:
     {
         auto meshB = std::dynamic_pointer_cast<SurfaceMesh>(collidingGeometryB);
         IMSTK_CHECK_FOR_VALID_GEOMETRIES(collidingGeometryA, meshB)
-            return std::make_shared<MeshToMeshBruteForceCD>(collidingGeometryA, meshB, colData);
+        return std::make_shared<MeshToMeshBruteForceCD>(collidingGeometryA, meshB, colData);
     }
     default:
     {
diff --git a/Source/CollisionDetection/imstkCDObjectFactory.h b/Source/CollisionDetection/imstkCDObjectFactory.h
index 63b48674454884ecf9b07b61df60c935d584997d..1f27972361706d8a9771a9a721057d8c7cabc54c 100644
--- a/Source/CollisionDetection/imstkCDObjectFactory.h
+++ b/Source/CollisionDetection/imstkCDObjectFactory.h
@@ -33,9 +33,8 @@ class Geometry;
 /// it will be added to an internal static octree for detecting collision
 /// \todo Other collision pair may be considered to use octree too
 ///
-extern std::shared_ptr<CollisionDetection>
-makeCollisionDetectionObject(const CollisionDetection::Type type,
-    std::shared_ptr<Geometry> collidingGeometryA,
-    std::shared_ptr<Geometry> collidingGeometryB,
-    std::shared_ptr<CollisionData> colData);
+extern std::shared_ptr<CollisionDetection> makeCollisionDetectionObject(const CollisionDetection::Type type,
+                                                                        std::shared_ptr<Geometry>      collidingGeometryA,
+                                                                        std::shared_ptr<Geometry>      collidingGeometryB,
+                                                                        std::shared_ptr<CollisionData> colData);
 }
diff --git a/Source/CollisionHandling/imstkCollisionHandling.cpp b/Source/CollisionHandling/imstkCollisionHandling.cpp
index d4aa1e92a38f11f003f78b9a91e3c675bafdb615..501b1b018518fca8faa8e4dc2f125ce3865b9274 100644
--- a/Source/CollisionHandling/imstkCollisionHandling.cpp
+++ b/Source/CollisionHandling/imstkCollisionHandling.cpp
@@ -29,6 +29,5 @@ CollisionHandling::CollisionHandling(const Type& type, const Side& side,
     m_type(type), m_side(side), m_colData(colData),
     m_computeNode(std::make_shared<ComputeNode>(std::bind(&CollisionHandling::processCollisionData, this), "CollisionHandling", true))
 {
-
 }
 }
\ No newline at end of file
diff --git a/Source/CollisionHandling/imstkPBDCollisionHandling.cpp b/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
index f7cc2cd602aacff3ad61d9b36bd9744364139638..31a2f57a37f9ffdf83f10da556ed648caf9bcfa6 100644
--- a/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
+++ b/Source/CollisionHandling/imstkPBDCollisionHandling.cpp
@@ -33,8 +33,8 @@ namespace imstk
 {
 PBDCollisionHandling::PBDCollisionHandling(const Side&                          side,
                                            const std::shared_ptr<CollisionData> colData,
-                                           std::shared_ptr<PbdObject> pbdObject1,
-                                           std::shared_ptr<PbdObject> pbdObject2) :
+                                           std::shared_ptr<PbdObject>           pbdObject1,
+                                           std::shared_ptr<PbdObject>           pbdObject2) :
     CollisionHandling(Type::PBD, side, colData),
     m_PbdObject1(pbdObject1),
     m_PbdObject2(pbdObject2),
diff --git a/Source/CollisionHandling/imstkPBDCollisionHandling.h b/Source/CollisionHandling/imstkPBDCollisionHandling.h
index 0f6fe3cd02fad783cfcee0cc84283decc5e2b1c1..3935e0db3328f3f4c93a77482f7b6fbae6a81cbe 100644
--- a/Source/CollisionHandling/imstkPBDCollisionHandling.h
+++ b/Source/CollisionHandling/imstkPBDCollisionHandling.h
@@ -49,8 +49,8 @@ public:
     ///
     PBDCollisionHandling(const Side&                          side,
                          const std::shared_ptr<CollisionData> colData,
-                         std::shared_ptr<PbdObject> pbdObject1,
-                         std::shared_ptr<PbdObject> pbdObject2);
+                         std::shared_ptr<PbdObject>           pbdObject1,
+                         std::shared_ptr<PbdObject>           pbdObject2);
 
     PBDCollisionHandling() = delete;
 
diff --git a/Source/CollisionHandling/imstkPenaltyCH.h b/Source/CollisionHandling/imstkPenaltyCH.h
index f6661bbd6f600e868a57b410c6ca6a00297eb9da..35b3724034bc90c0cb26bd141eb20640b372f121 100644
--- a/Source/CollisionHandling/imstkPenaltyCH.h
+++ b/Source/CollisionHandling/imstkPenaltyCH.h
@@ -79,7 +79,7 @@ public:
 private:
     std::shared_ptr<CollidingObject> m_object = nullptr; ///>
 
-    double m_stiffness = 5.0e5;                ///> Stiffness of contact
-    double m_damping   = 0.5;                  ///> Damping of the contact
+    double m_stiffness = 5.0e5;                          ///> Stiffness of contact
+    double m_damping   = 0.5;                            ///> Damping of the contact
 };
 }
diff --git a/Source/CollisionHandling/imstkSPHCollisionHandling.cpp b/Source/CollisionHandling/imstkSPHCollisionHandling.cpp
index 053043bf3e20802150175023bbb143993b250290..99670ce38a06b5cb74183a108335ebe09ccb8034 100644
--- a/Source/CollisionHandling/imstkSPHCollisionHandling.cpp
+++ b/Source/CollisionHandling/imstkSPHCollisionHandling.cpp
@@ -30,7 +30,7 @@ namespace imstk
 {
 SPHCollisionHandling::SPHCollisionHandling(const CollisionHandling::Side& side,
                                            std::shared_ptr<CollisionData> colData,
-                                           std::shared_ptr<SPHObject> obj) :
+                                           std::shared_ptr<SPHObject>     obj) :
     CollisionHandling(Type::SPH, side, colData), m_SPHObject(obj)
 {
     LOG_IF(FATAL, (!m_SPHObject)) << "Invalid SPH object";
@@ -54,9 +54,9 @@ SPHCollisionHandling::processCollisionData()
     ParallelUtils::parallelFor(m_colData->MAColData.getSize(),
         [&](const size_t idx)
         {
-            const MeshToAnalyticalCollisionDataElement& cd  = m_colData->MAColData[idx];
-            const uint32_t pidx = cd.nodeIdx;  // Fluid particle index
-            Vec3d n = cd.penetrationVector; // This vector should point into solid object
+            const MeshToAnalyticalCollisionDataElement& cd = m_colData->MAColData[idx];
+            const uint32_t pidx = cd.nodeIdx; // Fluid particle index
+            Vec3d n = cd.penetrationVector;   // This vector should point into solid object
 
             // Correct particle position
             state.getPositions()[pidx] -= n;
diff --git a/Source/CollisionHandling/imstkSPHCollisionHandling.h b/Source/CollisionHandling/imstkSPHCollisionHandling.h
index c7adb4a5ac5584bb48b08ac930d5eb108e35545d..e024c18d70375230690c0267cccf1c05a0c8e120 100644
--- a/Source/CollisionHandling/imstkSPHCollisionHandling.h
+++ b/Source/CollisionHandling/imstkSPHCollisionHandling.h
@@ -30,9 +30,9 @@ class SPHObject;
 class SPHCollisionHandling : public CollisionHandling
 {
 public:
-    SPHCollisionHandling(const Side& side,
+    SPHCollisionHandling(const Side&                    side,
                          std::shared_ptr<CollisionData> colData,
-                         std::shared_ptr<SPHObject> obj);
+                         std::shared_ptr<SPHObject>     obj);
 
     SPHCollisionHandling() = delete;
 
diff --git a/Source/Common/ComputeGraph/imstkComputeGraph.cpp b/Source/Common/ComputeGraph/imstkComputeGraph.cpp
index 8957da00501ccc1401be678a17830ac7ed4e2816..86fd4d2f523ca0850914cf9162de64f0ce2bd9ba 100644
--- a/Source/Common/ComputeGraph/imstkComputeGraph.cpp
+++ b/Source/Common/ComputeGraph/imstkComputeGraph.cpp
@@ -175,7 +175,7 @@ ComputeGraph::removeNode(std::shared_ptr<ComputeNode> node)
     }
 
     // Erase edges
-    ComputeNodeSet inputs = m_invAdjList[node];
+    ComputeNodeSet inputs  = m_invAdjList[node];
     ComputeNodeSet outputs = m_adjList[node];
     for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
     {
@@ -205,7 +205,7 @@ ComputeGraph::removeNodeAndFix(std::shared_ptr<ComputeNode> node)
     }
 
     // Erase edges
-    ComputeNodeSet inputs = m_invAdjList[node];
+    ComputeNodeSet inputs  = m_invAdjList[node];
     ComputeNodeSet outputs = m_adjList[node];
     for (ComputeNodeSet::iterator i = inputs.begin(); i != inputs.end(); i++)
     {
@@ -386,7 +386,7 @@ ComputeGraph::topologicalSort(std::shared_ptr<ComputeGraph> graph)
     // Compute the number of inputs to each node (we will remove these as we go)
     std::unordered_map<std::shared_ptr<ComputeNode>, size_t> numInputs;
 
-    const ComputeNodeAdjList& adjList = graph->getAdjList();
+    const ComputeNodeAdjList& adjList    = graph->getAdjList();
     const ComputeNodeAdjList& invAdjList = graph->getInvAdjList();
 
     for (ComputeNodeAdjList::const_iterator i = invAdjList.begin(); i != invAdjList.end(); i++)
@@ -449,7 +449,7 @@ ComputeGraph::resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph)
     std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
 
     const ComputeNodeAdjList& adjList = graph->getAdjList();
-    const ComputeNodeVector& nodes = graph->getNodes();
+    const ComputeNodeVector&  nodes   = graph->getNodes();
 
     // Compute the levels of each node via DFS
     std::unordered_map<std::shared_ptr<ComputeNode>, int> depths;
@@ -462,8 +462,8 @@ ComputeGraph::resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph)
         nodeStack.push(graph->getSource());
         while (!nodeStack.empty())
         {
-            std::shared_ptr<ComputeNode> currNode = nodeStack.top();
-            int currLevel = depths[currNode];
+            std::shared_ptr<ComputeNode> currNode  = nodeStack.top();
+            int                          currLevel = depths[currNode];
             nodeStack.pop();
 
             // Add children to stack if not yet visited
@@ -489,7 +489,9 @@ ComputeGraph::resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph)
     for (size_t i = 0; i < nodes.size(); i++)
     {
         if (nodes[i]->m_critical)
+        {
             critNodes.push_back(nodes[i]);
+        }
     }
 
     // Compute the critical adjacency list
@@ -538,8 +540,8 @@ ComputeGraph::resolveCriticalNodes(std::shared_ptr<ComputeGraph> graph)
         {
             std::shared_ptr<ComputeNode> destNode = critNodes[j];
             // If the edge doesn't exist, either way
-            if ((critAdjList.count(srcNode) == 0 || critAdjList.at(srcNode).count(destNode) == 0) &&
-                (critAdjList.count(destNode) == 0 || critAdjList.at(destNode).count(srcNode) == 0))
+            if ((critAdjList.count(srcNode) == 0 || critAdjList.at(srcNode).count(destNode) == 0)
+                && (critAdjList.count(destNode) == 0 || critAdjList.at(destNode).count(srcNode) == 0))
             {
                 // Add an edge between the critical nodes in the direction of levels
                 int leveli = depths[srcNode];
@@ -565,7 +567,9 @@ ComputeGraph::transitiveReduce(std::shared_ptr<ComputeGraph> graph)
 {
     // It's a bad idea to do this method if the graph is cyclic
     if (isCyclic(graph))
+    {
         return nullptr;
+    }
 
     std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
 
@@ -575,8 +579,8 @@ ComputeGraph::transitiveReduce(std::shared_ptr<ComputeGraph> graph)
     // For every edge in the graph
     for (ComputeNodeAdjList::const_iterator i = adjList.begin(); i != adjList.end(); i++)
     {
-        std::shared_ptr<ComputeNode> inputNode = i->first;
-        const ComputeNodeSet& outputNodes = i->second;
+        std::shared_ptr<ComputeNode> inputNode   = i->first;
+        const ComputeNodeSet&        outputNodes = i->second;
         for (ComputeNodeSet::const_iterator j = outputNodes.begin(); j != outputNodes.end(); j++)
         {
             std::shared_ptr<ComputeNode> outputNode = *j;
@@ -601,9 +605,9 @@ ComputeGraph::nonFunctionalPrune(std::shared_ptr<ComputeGraph> graph)
 {
     std::shared_ptr<ComputeGraph> results = std::make_shared<ComputeGraph>(*graph);
 
-    const ComputeNodeAdjList& adjList = results->getAdjList();
+    const ComputeNodeAdjList& adjList    = results->getAdjList();
     const ComputeNodeAdjList& invAdjList = results->getInvAdjList();
-    ComputeNodeVector& nodes = results->getNodes();
+    ComputeNodeVector&        nodes      = results->getNodes();
 
     for (size_t i = 0; i < nodes.size(); i++)
     {
@@ -619,7 +623,7 @@ ComputeGraph::nonFunctionalPrune(std::shared_ptr<ComputeGraph> graph)
         if (!node->isFunctional())
         {
             // Get copies of the inputs and outputs of the node
-            ComputeNodeSet inputs = invAdjList.at(node);
+            ComputeNodeSet inputs  = invAdjList.at(node);
             ComputeNodeSet outputs = adjList.at(node);
             if (inputs.size() == 1 && outputs.size() == 1)
             {
@@ -658,7 +662,7 @@ ComputeGraph::isCyclic(std::shared_ptr<ComputeGraph> graph)
 {
     // Brute force, DFS every node to find it again
     const ComputeNodeAdjList& adjList = graph->getAdjList();
-    const ComputeNodeVector& nodes = graph->getNodes();
+    const ComputeNodeVector&  nodes   = graph->getNodes();
     for (size_t i = 0; i < nodes.size(); i++)
     {
         std::unordered_set<std::shared_ptr<ComputeNode>> visitedNodes;
@@ -715,8 +719,8 @@ ComputeGraph::getUniqueNames(std::shared_ptr<ComputeGraph> graph, bool apply)
 {
     // Produce non colliding names
     std::unordered_map<std::shared_ptr<ComputeNode>, std::string> nodeNames;
-    std::unordered_map<std::string, int> names;
-    const ComputeNodeVector& nodes = graph->getNodes();
+    std::unordered_map<std::string, int>                          names;
+    const ComputeNodeVector&                                      nodes = graph->getNodes();
     for (size_t i = 0; i < nodes.size(); i++)
     {
         nodeNames[nodes[i]] = nodes[i]->m_name;
@@ -724,9 +728,9 @@ ComputeGraph::getUniqueNames(std::shared_ptr<ComputeGraph> graph, bool apply)
     }
     // Adjust names
     for (std::unordered_map<std::shared_ptr<ComputeNode>, std::string>::iterator it = nodeNames.begin();
-        it != nodeNames.end(); it++)
+         it != nodeNames.end(); it++)
     {
-        int nameIter = 0;
+        int         nameIter = 0;
         std::string currName = it->second;
         // If we can find a node with this name, increment name counter and try again
         while (names[currName] > 1)
@@ -794,7 +798,7 @@ ComputeGraph::getCriticalPath(std::shared_ptr<ComputeGraph> graph)
     std::unordered_map<std::shared_ptr<ComputeNode>, double> times = getTimes(graph);
 
     // Now backtrack to acquire the path of longest duration
-    const ComputeNodeAdjList& invAdjList = graph->getInvAdjList();
+    const ComputeNodeAdjList&               invAdjList = graph->getInvAdjList();
     std::list<std::shared_ptr<ComputeNode>> results;
     {
         std::shared_ptr<ComputeNode> currNode = graph->getSink();
@@ -803,7 +807,7 @@ ComputeGraph::getCriticalPath(std::shared_ptr<ComputeGraph> graph)
         {
             results.push_front(currNode);
             std::shared_ptr<ComputeNode> longestNode = nullptr;
-            double maxTime = 0.0;
+            double                       maxTime     = 0.0;
 
             // For every parent
             if (invAdjList.count(currNode) != 0)
@@ -814,7 +818,7 @@ ComputeGraph::getCriticalPath(std::shared_ptr<ComputeGraph> graph)
                     std::shared_ptr<ComputeNode> parentNode = *i;
                     if (times[parentNode] >= maxTime)
                     {
-                        maxTime = times[parentNode];
+                        maxTime     = times[parentNode];
                         longestNode = parentNode;
                     }
                 }
diff --git a/Source/Common/ComputeGraph/imstkComputeGraph.h b/Source/Common/ComputeGraph/imstkComputeGraph.h
index ae2422f16a821681c605afac7783040e014a6da4..eb249c012faaceb4900f31d81ea2420281c47300 100644
--- a/Source/Common/ComputeGraph/imstkComputeGraph.h
+++ b/Source/Common/ComputeGraph/imstkComputeGraph.h
@@ -66,7 +66,7 @@ public:
     ///
     ComputeNodeAdjList& getInvAdjList() { return m_invAdjList; }
 
-    // Node operations
+// Node operations
 public:
     ///
     /// \brief Linear search for node by name within this graph
@@ -130,7 +130,7 @@ public:
     ///
     void insertBefore(std::shared_ptr<ComputeNode> refNode, std::shared_ptr<ComputeNode> newNode);
 
-    // Edge operations
+// Edge operations
 public:
     ///
     /// \brief Returns whether or not this graph contains the given directed edge
@@ -173,7 +173,7 @@ public:
         m_invAdjList.clear();
     }
 
-    // Graph algorithms, todo: Move into filtering module
+// Graph algorithms, todo: Move into filtering module
 public:
     ///
     /// \brief Graph sum, shared references are considered identical nodes, source/sink of results invalidated/nullptr
@@ -232,7 +232,7 @@ public:
 
 protected:
     ComputeNodeVector  m_nodes;
-    ComputeNodeAdjList m_adjList; ///> This gives the outputs of every node
+    ComputeNodeAdjList m_adjList;    ///> This gives the outputs of every node
     ComputeNodeAdjList m_invAdjList; ///> This gives the inputs of every node
     // inv adjlist?
 
diff --git a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp
index 895416cda7d96344ef68cbacf830c213dbd1b903..6159308842faaec9bf1a811ff2f3f1c0e644504a 100644
--- a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp
+++ b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.cpp
@@ -64,11 +64,11 @@ ComputeGraphVizWriter::write()
         critPath = ComputeGraph::getCriticalPath(m_inputGraph);
     }
     auto edgeExists = [&](const std::shared_ptr<ComputeNode>& a, const std::shared_ptr<ComputeNode>& b)
-    {
-        ComputeNodeList::iterator srcNode = std::find(critPath.begin(), critPath.end(), a);
-        // If srcNode was found and the next node is b
-        return (srcNode != critPath.end() && *std::next(srcNode) == b);
-    };
+                      {
+                          ComputeNodeList::iterator srcNode = std::find(critPath.begin(), critPath.end(), a);
+                          // If srcNode was found and the next node is b
+                          return (srcNode != critPath.end() && *std::next(srcNode) == b);
+                      };
 
     // Compute completion times of each node
     if (m_writeTimes)
diff --git a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h
index 8d07c9692662427667ccbe7973a8c28e0002d32c..3d42da8188b5c0ad0a1af5b4a1ad19b3ef5ac9fc 100644
--- a/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h
+++ b/Source/Common/ComputeGraph/imstkComputeGraphVizWriter.h
@@ -75,7 +75,7 @@ public:
 
 private:
     std::shared_ptr<ComputeGraph> m_inputGraph = nullptr;
-    std::string m_fileName = "";
+    std::string m_fileName       = "";
     bool m_highlightCriticalPath = false;
     bool m_writeTimes = false;
 };
diff --git a/Source/Common/ComputeGraph/imstkComputeNode.h b/Source/Common/ComputeGraph/imstkComputeNode.h
index 31c519f869a050ab9f36d971b953df1539f5b50c..daa570e1675484f0512b7f071bfd7e9076b04b39 100644
--- a/Source/Common/ComputeGraph/imstkComputeNode.h
+++ b/Source/Common/ComputeGraph/imstkComputeNode.h
@@ -55,11 +55,11 @@ public:
     virtual void execute();
 
 public:
-    std::string m_name = "none";
-    bool m_enabled     = true;
-    bool m_critical    = false;
+    std::string m_name   = "none";
+    bool   m_enabled     = true;
+    bool   m_critical    = false;
     double m_elapsedTime = 0.0;
-    bool m_enableBenchmarking = false;
+    bool   m_enableBenchmarking = false;
 
 protected:
     std::function<void()> m_func = nullptr; ///> Don't allow user to call directly (must use execute)
diff --git a/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h
index c8e1f2156b32f3787fa1f8b4b9d31aae19062da6..ba89cb3fbfb36832c1865861655f6ba0fd4d8ce5 100644
--- a/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h
+++ b/Source/Common/ComputeGraph/imstkSequentialComputeGraphController.h
@@ -45,6 +45,6 @@ public:
 
 private:
     // The current nodes to execute, ordered
-    std::shared_ptr<std::list<std::shared_ptr<ComputeNode>>>  m_executionOrderedNodes;
+    std::shared_ptr<std::list<std::shared_ptr<ComputeNode>>> m_executionOrderedNodes;
 };
 };
\ No newline at end of file
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index 2a29659ae7f87a1dbd6a5c9adf87398d250a50a4..82be1297d64a5ebe0b6b5e809a9f371ae3bbedfa 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -32,7 +32,6 @@
 #include "imstkIsotropicHyperelasticFEMForceModel.h"
 #include "imstkLinearFEMForceModel.h"
 #include "imstkMath.h"
-#include "imstkNewtonSolver.h"
 #include "imstkSolverBase.h"
 #include "imstkStVKForceModel.h"
 #include "imstkTimeIntegrator.h"
@@ -60,7 +59,7 @@ FEMDeformableBodyModel::FEMDeformableBodyModel() :
         Geometry::Type::HexahedralMesh
     };
 
-    m_solveNode      = addFunction("FEMModel_Solve", [&]() { getSolver()->solve(); });
+    m_solveNode = addFunction("FEMModel_Solve", [&]() { getSolver()->solve(); });
 }
 
 FEMDeformableBodyModel::~FEMDeformableBodyModel()
@@ -167,7 +166,7 @@ FEMDeformableBodyModel::initialize()
     if (m_solver == nullptr)
     {
         // Create a nonlinear system
-        auto nlSystem = std::make_shared<NonLinearSystem>(getFunction(), getFunctionGradient());
+        auto nlSystem = std::make_shared<NonLinearSystem<SparseMatrixd>>(getFunction(), getFunctionGradient());
 
         nlSystem->setUnknownVector(getUnknownVec());
         nlSystem->setUpdateFunction(getUpdateFunction());
@@ -183,7 +182,7 @@ FEMDeformableBodyModel::initialize()
         }
 
         // Create a non-linear solver and add to the scene
-        auto nlSolver = std::make_shared<NewtonSolver>();
+        auto nlSolver = std::make_shared<NewtonSolver<SparseMatrixd>>();
         nlSolver->setLinearSolver(linSolver);
         nlSolver->setSystem(nlSystem);
         setSolver(nlSolver);
@@ -289,7 +288,7 @@ FEMDeformableBodyModel::initializeForceModel()
     const double g = m_FEModelConfig->m_gravity;
     // Since vega 4.0 doesn't add gravity correcntly in all cases, we do it ourselves; see \ref initializeGravityForce
     // const bool   isGravityPresent = (g > 0) ? true : false;
-    const bool   isGravityPresent = false;
+    const bool isGravityPresent = false;
 
     m_numDOF = (size_t)m_vegaPhysicsMesh->getNumVertices() * 3;
 
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index 3eeafd96e8dff4c2b55a59004bb634026a977a0f..53209f3421c17b91a6bac7ae607aff9dbdc11f0d 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -25,6 +25,7 @@
 #include "imstkDynamicalModel.h"
 #include "imstkInternalForceModelTypes.h"
 #include "imstkNonLinearSystem.h"
+#include "imstkNewtonSolver.h"
 
 // vega
 #include "sparseMatrix.h"
@@ -37,7 +38,6 @@ class VolumetricMesh;
 namespace imstk
 {
 class InternalForceModel;
-class NewtonSolver;
 class TimeIntegrator;
 class SolverBase;
 class VegaMeshIO;
@@ -279,9 +279,9 @@ protected:
 
 protected:
     std::shared_ptr<SolverBase> m_solver = nullptr;
-    std::shared_ptr<InternalForceModel> m_internalForceModel = nullptr;       ///> Mathematical model for intenal forces
-    std::shared_ptr<TimeIntegrator>     m_timeIntegrator     = nullptr;       ///> Time integrator
-    std::shared_ptr<NonLinearSystem>    m_nonLinearSystem    = nullptr;       ///> Nonlinear system resulting from TI and force model
+    std::shared_ptr<InternalForceModel> m_internalForceModel = nullptr;          ///> Mathematical model for intenal forces
+    std::shared_ptr<TimeIntegrator>     m_timeIntegrator     = nullptr;          ///> Time integrator
+    std::shared_ptr<NonLinearSystem<SparseMatrixd>> m_nonLinearSystem = nullptr; ///> Nonlinear system resulting from TI and force model
 
     std::shared_ptr<FEMModelConfig> m_FEModelConfig = nullptr;
 
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index c03caa4e302e36393ef67acb0db581143303448e..9f8293b68dd86c3bf6f4ede5f5330aec06c0f6bd 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -57,10 +57,10 @@ PbdModel::PbdModel() : DynamicalModel(DynamicalModelType::PositionBasedDynamics)
     };
 
     // Setup PBD compute nodes
-    m_integrationPositionNode = addFunction("PbdModel_IntegratePosition", std::bind(&PbdModel::integratePosition, this));
+    m_integrationPositionNode     = addFunction("PbdModel_IntegratePosition", std::bind(&PbdModel::integratePosition, this));
     m_updateCollisionGeometryNode = addFunction("PbdModel_UpdateCollisionGeometry", std::bind(&PbdModel::updatePhysicsGeometry, this));
-    m_solveConstraintsNode    = addFunction("PbdModel_SolveConstraints", [&]() { m_pbdSolver->solve(); }); // Avoids rebinding on solver swap
-    m_updateVelocityNode      = addFunction("PbdModel_UpdateVelocity", std::bind(&PbdModel::updateVelocity, this));
+    m_solveConstraintsNode = addFunction("PbdModel_SolveConstraints", [&]() { m_pbdSolver->solve(); });    // Avoids rebinding on solver swap
+    m_updateVelocityNode   = addFunction("PbdModel_UpdateVelocity", std::bind(&PbdModel::updateVelocity, this));
 }
 
 void
diff --git a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
index 7f34a3c9250ac4c974e381b6bf6ebac1a076f9d7..c21e2f04d92e17e753645e91461618ee033d671b 100644
--- a/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -276,9 +276,9 @@ protected:
 
 protected:
     // Computational Nodes
-    std::shared_ptr<ComputeNode> m_integrationPositionNode = nullptr;
+    std::shared_ptr<ComputeNode> m_integrationPositionNode     = nullptr;
     std::shared_ptr<ComputeNode> m_updateCollisionGeometryNode = nullptr;
-    std::shared_ptr<ComputeNode> m_solveConstraintsNode    = nullptr;
-    std::shared_ptr<ComputeNode> m_updateVelocityNode      = nullptr;
+    std::shared_ptr<ComputeNode> m_solveConstraintsNode = nullptr;
+    std::shared_ptr<ComputeNode> m_updateVelocityNode   = nullptr;
 };
 } // imstk
diff --git a/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp b/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp
index 66cd3049befee8536c02adcb3f1c07bb3d78e800..f155fca619bff742b00f3660e72c5dc8100e66ea 100644
--- a/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp
+++ b/Source/Rendering/VTKRenderer/imstkVTKRenderer.cpp
@@ -415,7 +415,6 @@ VTKRenderer::getAxesVisibility() const
     return m_AxesActor->GetVisibility();
 }
 
-
 void
 VTKRenderer::setBenchmarkTable(const std::unordered_map<std::string, double>& nameToElapsedTimes)
 {
@@ -425,7 +424,7 @@ VTKRenderer::setBenchmarkTable(const std::unordered_map<std::string, double>& na
         [](const std::pair<std::string, double>& a, const std::pair<std::string, double>& b) { return a.second < b.second; });
 
     // Construct vtkTable from provided data
-    vtkSmartPointer<vtkDoubleArray> xIndices = vtkDoubleArray::SafeDownCast(m_benchmarkTable->GetColumn(0));
+    vtkSmartPointer<vtkDoubleArray> xIndices      = vtkDoubleArray::SafeDownCast(m_benchmarkTable->GetColumn(0));
     vtkSmartPointer<vtkDoubleArray> yElapsedTimes = vtkDoubleArray::SafeDownCast(m_benchmarkTable->GetColumn(1));
     vtkSmartPointer<vtkStringArray> labels = vtkStringArray::SafeDownCast(m_benchmarkTable->GetColumn(2));
 
@@ -487,7 +486,6 @@ VTKRenderer::getBenchmarkTableVisibility() const
     return m_benchmarkChartActor->GetVisibility();
 }
 
-
 void
 VTKRenderer::updateSceneCamera(std::shared_ptr<Camera> imstkCam)
 {
diff --git a/Source/Rendering/VTKRenderer/imstkVTKRenderer.h b/Source/Rendering/VTKRenderer/imstkVTKRenderer.h
index e55da7653fa3e7de6953b02981dc7a70270d0923..b7e378502953c24d1955a84dd07b63b54d9ce1d2 100644
--- a/Source/Rendering/VTKRenderer/imstkVTKRenderer.h
+++ b/Source/Rendering/VTKRenderer/imstkVTKRenderer.h
@@ -168,7 +168,7 @@ protected:
     std::vector<vtkOpenVRCameraPose> m_camPos;
 #endif
 
-    vtkSmartPointer<vtkChartXY> m_benchmarkChart = nullptr;
+    vtkSmartPointer<vtkChartXY>      m_benchmarkChart      = nullptr;
     vtkSmartPointer<vtkContextActor> m_benchmarkChartActor = nullptr;
     vtkSmartPointer<vtkTable> m_benchmarkTable = nullptr;
     vtkPlotBar* m_benchmarkPlot = nullptr;
diff --git a/Source/Scene/imstkCollisionGraph.cpp b/Source/Scene/imstkCollisionGraph.cpp
index afee115674c523cfcb5e6227a392a9cdde547496..6805e2faa8b2f64aa9d9dab2681bdf7f4978a515 100644
--- a/Source/Scene/imstkCollisionGraph.cpp
+++ b/Source/Scene/imstkCollisionGraph.cpp
@@ -154,9 +154,9 @@ CollisionGraph::removeInteractionPair(SceneObjectPtr A, SceneObjectPtr B)
 bool
 CollisionGraph::removeInteractionPair(ObjectInteractionPtr intPair)
 {
-   /* std::shared_ptr<CollidingObject> obj1 = std::static_pointer_cast<CollidingObject>(intPair->getObjectsPair().first);
-    std::shared_ptr<CollidingObject> obj2 = std::static_pointer_cast<CollidingObject>(intPair->getObjectsPair().second);
-    return this->removeInteractionPair(obj1, obj2);*/
+    /* std::shared_ptr<CollidingObject> obj1 = std::static_pointer_cast<CollidingObject>(intPair->getObjectsPair().first);
+     std::shared_ptr<CollidingObject> obj2 = std::static_pointer_cast<CollidingObject>(intPair->getObjectsPair().second);
+     return this->removeInteractionPair(obj1, obj2);*/
     return true;
 }
 
diff --git a/Source/Scene/imstkCollisionGraph.h b/Source/Scene/imstkCollisionGraph.h
index 9a723c994dc84835ede2a467fdc761d1bbf654e6..0438862499b7981b65f562ac18073ff55a04ba0e 100644
--- a/Source/Scene/imstkCollisionGraph.h
+++ b/Source/Scene/imstkCollisionGraph.h
@@ -27,7 +27,6 @@
 
 namespace imstk
 {
-
 class SceneObject;
 class ObjectInteractionPair;
 
@@ -39,8 +38,8 @@ class ObjectInteractionPair;
 class CollisionGraph
 {
 public:
-    using SceneObjectPtr        = std::shared_ptr<SceneObject>;
-    using ObjectInteractionPtr  = std::shared_ptr<ObjectInteractionPair>;
+    using SceneObjectPtr       = std::shared_ptr<SceneObject>;
+    using ObjectInteractionPtr = std::shared_ptr<ObjectInteractionPair>;
 
     ///
     /// \brief Default constructor
@@ -93,7 +92,7 @@ public:
     const std::unordered_map<SceneObjectPtr, std::vector<ObjectInteractionPtr>>& getInteractionPairMap() const;
 
 protected:
-    std::vector<ObjectInteractionPtr> m_interactionPairs;                                           ///< All interaction pairs in the collision graph
+    std::vector<ObjectInteractionPtr> m_interactionPairs;                                       ///< All interaction pairs in the collision graph
     std::unordered_map<SceneObjectPtr, std::vector<ObjectInteractionPtr>> m_interactionPairMap; ///< Map of interaction pairs per colliding object
 };
 }
diff --git a/Source/Scene/imstkCollisionPair.cpp b/Source/Scene/imstkCollisionPair.cpp
index ebb88a1b2a65c14ebf4db41e6d5c1f41c296b978..fe774d9b44dba870648372c9edf67aca1523a88d 100644
--- a/Source/Scene/imstkCollisionPair.cpp
+++ b/Source/Scene/imstkCollisionPair.cpp
@@ -32,159 +32,159 @@ limitations under the License.
 namespace imstk
 {
 CollisionPair::CollisionPair(std::shared_ptr<CollidingObject> objA,
-	std::shared_ptr<CollidingObject> objB) : ObjectInteractionPair(objA, objB)
+                             std::shared_ptr<CollidingObject> objB) : ObjectInteractionPair(objA, objB)
 {
 }
 
 CollisionPair::CollisionPair(std::shared_ptr<CollidingObject>    objA,
-	std::shared_ptr<CollidingObject>    objB,
-	std::shared_ptr<CollisionDetection> cd,
-	std::shared_ptr<CollisionHandling>  chA,
-	std::shared_ptr<CollisionHandling>  chB) : ObjectInteractionPair(objA, objB)
+                             std::shared_ptr<CollidingObject>    objB,
+                             std::shared_ptr<CollisionDetection> cd,
+                             std::shared_ptr<CollisionHandling>  chA,
+                             std::shared_ptr<CollisionHandling>  chB) : ObjectInteractionPair(objA, objB)
 {
-	setCollisionDetection(cd);
-
-	if (chA != nullptr)
-	{
-		setCollisionHandlingA(chA);
-	}
-	if (chB != nullptr)
-	{
-		setCollisionHandlingB(chB);
-	}
+    setCollisionDetection(cd);
+
+    if (chA != nullptr)
+    {
+        setCollisionHandlingA(chA);
+    }
+    if (chB != nullptr)
+    {
+        setCollisionHandlingB(chB);
+    }
 }
 
 CollisionPair::CollisionPair(std::shared_ptr<CollidingObject> objA, std::shared_ptr<CollidingObject> objB,
-	std::shared_ptr<CollisionDetection> cd,
-	std::shared_ptr<CollisionHandling> chAB) : ObjectInteractionPair(objA, objB)
+                             std::shared_ptr<CollisionDetection> cd,
+                             std::shared_ptr<CollisionHandling> chAB) : ObjectInteractionPair(objA, objB)
 {
-	setCollisionDetection(cd);
+    setCollisionDetection(cd);
 
-	if (chAB != nullptr)
-	{
-		setCollisionHandlingAB(chAB);
-	}
+    if (chAB != nullptr)
+    {
+        setCollisionHandlingAB(chAB);
+    }
 }
 
 void
 CollisionPair::setCollisionDetection(std::shared_ptr<CollisionDetection> colDetect)
 {
-	m_colDetect = colDetect;
-	m_collisionDetectionNode = m_interactionFunction = m_colDetect->getComputeNode();
-	m_collisionDetectionNode->m_name = getObjectsPair().first->getName() + "_" + getObjectsPair().second->getName() + "_CollisionDetection";
-	m_colData = m_colDetect->getCollisionData();
+    m_colDetect = colDetect;
+    m_collisionDetectionNode = m_interactionFunction = m_colDetect->getComputeNode();
+    m_collisionDetectionNode->m_name = getObjectsPair().first->getName() + "_" + getObjectsPair().second->getName() + "_CollisionDetection";
+    m_colData = m_colDetect->getCollisionData();
 }
 
 void
 CollisionPair::setCollisionHandlingA(std::shared_ptr<CollisionHandling> colHandlingA)
 {
-	m_colHandlingA = colHandlingA;
-	m_collisionHandleANode = m_colHandlingA->getComputeNode();
-	m_collisionHandleANode->m_name = getObjectsPair().first->getName() + "_CollisionHandling";
+    m_colHandlingA = colHandlingA;
+    m_collisionHandleANode = m_colHandlingA->getComputeNode();
+    m_collisionHandleANode->m_name = getObjectsPair().first->getName() + "_CollisionHandling";
 }
 
 void
 CollisionPair::setCollisionHandlingB(std::shared_ptr<CollisionHandling> colHandlingB)
 {
-	m_colHandlingB = colHandlingB;
-	m_collisionHandleBNode = m_colHandlingB->getComputeNode();
-	m_collisionHandleBNode->m_name = getObjectsPair().second->getName() + "_CollisionHandling";
+    m_colHandlingB = colHandlingB;
+    m_collisionHandleBNode = m_colHandlingB->getComputeNode();
+    m_collisionHandleBNode->m_name = getObjectsPair().second->getName() + "_CollisionHandling";
 }
 
 void
 CollisionPair::setCollisionHandlingAB(std::shared_ptr<CollisionHandling> colHandlingAB)
 {
-	m_colHandlingA = m_colHandlingB = colHandlingAB;
-	m_collisionHandleANode = m_collisionHandleBNode = colHandlingAB->getComputeNode();
-	m_collisionHandleANode->m_name = getObjectsPair().first->getName() + '_' + getObjectsPair().second->getName() + "_CollisionHandling";
+    m_colHandlingA = m_colHandlingB = colHandlingAB;
+    m_collisionHandleANode = m_collisionHandleBNode = colHandlingAB->getComputeNode();
+    m_collisionHandleANode->m_name = getObjectsPair().first->getName() + '_' + getObjectsPair().second->getName() + "_CollisionHandling";
 }
 
 void
 CollisionPair::modifyComputeGraph()
 {
-	std::shared_ptr<ComputeGraph> computeGraphA = m_objects.first->getComputeGraph();
-	std::shared_ptr<ComputeGraph> computeGraphB = m_objects.second->getComputeGraph();
-
-	// If nothing was added to the input/output list use default collision location
-	if ((m_computeNodeInputs.first.size() == 0) && (m_computeNodeInputs.second.size() == 0) &&
-		(m_computeNodeOutputs.first.size() == 0) && (m_computeNodeOutputs.second.size() == 0))
-	{
-		m_computeNodeInputs.first.clear();
-		m_computeNodeInputs.second.clear();
-		m_computeNodeOutputs.first.clear();
-		m_computeNodeOutputs.second.clear();
-
-		// Default location is the first node in the SceneObject
-		m_computeNodeInputs.first.push_back(computeGraphA->getSource());
-		m_computeNodeInputs.first.push_back(computeGraphB->getSource());
-
-		m_computeNodeInputs.first.push_back(m_objects.first->getUpdateNode());
-		m_computeNodeInputs.first.push_back(m_objects.second->getUpdateNode());
-	}
-
-	// Add all the nodes to the graph
-	computeGraphA->addNode(m_collisionDetectionNode);
-	computeGraphB->addNode(m_collisionDetectionNode);
-	if (m_collisionHandleANode != nullptr)
-	{
-		computeGraphA->addNode(m_collisionHandleANode);
-	}
-	if (m_collisionHandleBNode != nullptr)
-	{
-		computeGraphB->addNode(m_collisionHandleBNode);
-	}
-
-	// Add the edges
-	{
-		// Connect inputA's->CD
-		for (size_t i = 0; i < m_computeNodeInputs.first.size(); i++)
-		{
-			computeGraphA->addEdge(m_computeNodeInputs.first[i], m_collisionDetectionNode);
-		}
-
-		// Connect inputB's->CD
-		for (size_t i = 0; i < m_computeNodeInputs.second.size(); i++)
-		{
-			computeGraphA->addEdge(m_computeNodeInputs.second[i], m_collisionDetectionNode);
-		}
-	}
-
-	// Now connect CD to CHA/CHB/CHAB
-	// This also works for the case CHA = CHB = CHAB
-	{
-		if (m_collisionHandleANode != nullptr)
-		{
-			computeGraphA->addEdge(m_collisionDetectionNode, m_collisionHandleANode);
-		}
-		if (m_collisionHandleBNode != nullptr)
-		{
-			computeGraphB->addEdge(m_collisionDetectionNode, m_collisionHandleBNode);
-		}
-	}
-
-	// Connect either CD or CHA/CHAB to outputA's
-	for (size_t i = 0; i < m_computeNodeOutputs.first.size(); i++)
-	{
-		if (m_collisionHandleANode != nullptr)
-		{
-			computeGraphA->addEdge(m_collisionHandleANode, m_computeNodeOutputs.first[i]);
-		}
-		else
-		{
-			computeGraphA->addEdge(m_collisionDetectionNode, m_computeNodeOutputs.first[i]);
-		}
-	}
-	// Connect eitehr CD or CHB/CHAB to outputB's
-	for (size_t i = 0; i < m_computeNodeOutputs.second.size(); i++)
-	{
-		if (m_collisionHandleBNode != nullptr)
-		{
-			computeGraphB->addEdge(m_collisionHandleBNode, m_computeNodeOutputs.second[i]);
-		}
-		else
-		{
-			computeGraphB->addEdge(m_collisionDetectionNode, m_computeNodeOutputs.second[i]);
-		}
-	}
+    std::shared_ptr<ComputeGraph> computeGraphA = m_objects.first->getComputeGraph();
+    std::shared_ptr<ComputeGraph> computeGraphB = m_objects.second->getComputeGraph();
+
+    // If nothing was added to the input/output list use default collision location
+    if ((m_computeNodeInputs.first.size() == 0) && (m_computeNodeInputs.second.size() == 0)
+        && (m_computeNodeOutputs.first.size() == 0) && (m_computeNodeOutputs.second.size() == 0))
+    {
+        m_computeNodeInputs.first.clear();
+        m_computeNodeInputs.second.clear();
+        m_computeNodeOutputs.first.clear();
+        m_computeNodeOutputs.second.clear();
+
+        // Default location is the first node in the SceneObject
+        m_computeNodeInputs.first.push_back(computeGraphA->getSource());
+        m_computeNodeInputs.first.push_back(computeGraphB->getSource());
+
+        m_computeNodeInputs.first.push_back(m_objects.first->getUpdateNode());
+        m_computeNodeInputs.first.push_back(m_objects.second->getUpdateNode());
+    }
+
+    // Add all the nodes to the graph
+    computeGraphA->addNode(m_collisionDetectionNode);
+    computeGraphB->addNode(m_collisionDetectionNode);
+    if (m_collisionHandleANode != nullptr)
+    {
+        computeGraphA->addNode(m_collisionHandleANode);
+    }
+    if (m_collisionHandleBNode != nullptr)
+    {
+        computeGraphB->addNode(m_collisionHandleBNode);
+    }
+
+    // Add the edges
+    {
+        // Connect inputA's->CD
+        for (size_t i = 0; i < m_computeNodeInputs.first.size(); i++)
+        {
+            computeGraphA->addEdge(m_computeNodeInputs.first[i], m_collisionDetectionNode);
+        }
+
+        // Connect inputB's->CD
+        for (size_t i = 0; i < m_computeNodeInputs.second.size(); i++)
+        {
+            computeGraphA->addEdge(m_computeNodeInputs.second[i], m_collisionDetectionNode);
+        }
+    }
+
+    // Now connect CD to CHA/CHB/CHAB
+    // This also works for the case CHA = CHB = CHAB
+    {
+        if (m_collisionHandleANode != nullptr)
+        {
+            computeGraphA->addEdge(m_collisionDetectionNode, m_collisionHandleANode);
+        }
+        if (m_collisionHandleBNode != nullptr)
+        {
+            computeGraphB->addEdge(m_collisionDetectionNode, m_collisionHandleBNode);
+        }
+    }
+
+    // Connect either CD or CHA/CHAB to outputA's
+    for (size_t i = 0; i < m_computeNodeOutputs.first.size(); i++)
+    {
+        if (m_collisionHandleANode != nullptr)
+        {
+            computeGraphA->addEdge(m_collisionHandleANode, m_computeNodeOutputs.first[i]);
+        }
+        else
+        {
+            computeGraphA->addEdge(m_collisionDetectionNode, m_computeNodeOutputs.first[i]);
+        }
+    }
+    // Connect eitehr CD or CHB/CHAB to outputB's
+    for (size_t i = 0; i < m_computeNodeOutputs.second.size(); i++)
+    {
+        if (m_collisionHandleBNode != nullptr)
+        {
+            computeGraphB->addEdge(m_collisionHandleBNode, m_computeNodeOutputs.second[i]);
+        }
+        else
+        {
+            computeGraphB->addEdge(m_collisionDetectionNode, m_computeNodeOutputs.second[i]);
+        }
+    }
 }
 }
diff --git a/Source/Scene/imstkInteractionPair.h b/Source/Scene/imstkInteractionPair.h
index 9db4f0ad8acab9fea6b59b40980869c2ce3c0a95..4997ff7839e418204266bb5a3c990c83a412c52b 100644
--- a/Source/Scene/imstkInteractionPair.h
+++ b/Source/Scene/imstkInteractionPair.h
@@ -36,7 +36,7 @@ class ComputeNode;
 class InteractionPair
 {
 public:
-    using Inputs = std::pair<std::vector<std::shared_ptr<ComputeNode>>, std::vector<std::shared_ptr<ComputeNode>>>;
+    using Inputs  = std::pair<std::vector<std::shared_ptr<ComputeNode>>, std::vector<std::shared_ptr<ComputeNode>>>;
     using Outputs = std::pair<std::vector<std::shared_ptr<ComputeNode>>, std::vector<std::shared_ptr<ComputeNode>>>;
 
 public:
@@ -48,7 +48,7 @@ public:
     const Outputs& getComputeNodeOutputs() const { return m_computeNodeOutputs; }
 
 protected:
-    Inputs m_computeNodeInputs;                         ///> The interacting nodes
+    Inputs  m_computeNodeInputs;                        ///> The interacting nodes
     Outputs m_computeNodeOutputs;                       ///> The interacting nodes
     std::shared_ptr<ComputeNode> m_interactionFunction; ///> Function to execute on interaction
 };
diff --git a/Source/Scene/imstkObjectInteractionFactory.h b/Source/Scene/imstkObjectInteractionFactory.h
index 9783b238bfe9d660d83337c71b4d845165fd0685..65f25ccc6a8657670803dfb92bbfc37f320f014d 100644
--- a/Source/Scene/imstkObjectInteractionFactory.h
+++ b/Source/Scene/imstkObjectInteractionFactory.h
@@ -61,15 +61,16 @@ enum class InteractionType
 /// \brief Factory for InteractionPairs
 ///
 template<typename ObjectType1, typename ObjectType2>
-extern std::shared_ptr<ObjectInteractionPair> makeObjectInteractionPair(
+extern std::shared_ptr<ObjectInteractionPair>
+makeObjectInteractionPair(
     std::shared_ptr<ObjectType1> obj1, std::shared_ptr<ObjectType2> obj2,
     InteractionType intType, CollisionDetection::Type cdType)
 {
     std::shared_ptr<ObjectInteractionPair> results = nullptr;
     if (intType == InteractionType::PbdObjToPbdObj_Collision)
     {
-        if (std::is_base_of<PbdObject, ObjectType1>::value &&
-            std::is_base_of<PbdObject, ObjectType2>::value)
+        if (std::is_base_of<PbdObject, ObjectType1>::value
+            && std::is_base_of<PbdObject, ObjectType2>::value)
         {
             results = std::make_shared<PbdObjectCollisionPair>(
                 std::dynamic_pointer_cast<PbdObject>(obj1),
@@ -88,27 +89,27 @@ extern std::shared_ptr<ObjectInteractionPair> makeObjectInteractionPair(
     }*/
     else if (intType == InteractionType::SphObjToCollidingObj_Collision)
     {
-        if (std::is_base_of<SPHObject, ObjectType1>::value &&
-            std::is_base_of<CollidingObject, ObjectType2>::value)
+        if (std::is_base_of<SPHObject, ObjectType1>::value
+            && std::is_base_of<CollidingObject, ObjectType2>::value)
         {
             // Setup CD, and collision data
-            std::shared_ptr<CollisionData> colData = std::make_shared<CollisionData>();
+            std::shared_ptr<CollisionData>      colData   = std::make_shared<CollisionData>();
             std::shared_ptr<CollisionDetection> colDetect = makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), colData);
 
             // Setup the handler
             std::shared_ptr<SPHCollisionHandling> colHandler =
                 std::make_shared<SPHCollisionHandling>(CollisionHandling::Side::A, colData, std::dynamic_pointer_cast<SPHObject>(obj1));
-           
+
             results = std::make_shared<CollisionPair>(obj1, obj2, colDetect, colHandler, nullptr);
         }
     }
     else if (intType == InteractionType::FemObjToCollidingObj_NodalPicking)
     {
-        if (std::is_base_of<FeDeformableObject, ObjectType1>::value &&
-            std::is_base_of<CollidingObject, ObjectType2>::value)
+        if (std::is_base_of<FeDeformableObject, ObjectType1>::value
+            && std::is_base_of<CollidingObject, ObjectType2>::value)
         {
             // Setup CD, and collision data
-            std::shared_ptr<CollisionData> colData = std::make_shared<CollisionData>();
+            std::shared_ptr<CollisionData>      colData   = std::make_shared<CollisionData>();
             std::shared_ptr<CollisionDetection> colDetect = makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), colData);
 
             // Setup the handler
@@ -120,11 +121,11 @@ extern std::shared_ptr<ObjectInteractionPair> makeObjectInteractionPair(
     }
     else if (intType == InteractionType::FemObjToCollidingObj_PenaltyForce)
     {
-        if ((std::is_base_of<CollidingObject, ObjectType1>::value) &&
-            std::is_base_of<CollidingObject, ObjectType2>::value)
+        if ((std::is_base_of<CollidingObject, ObjectType1>::value)
+            && std::is_base_of<CollidingObject, ObjectType2>::value)
         {
             // Setup CD, and collision data
-            std::shared_ptr<CollisionData> colData = std::make_shared<CollisionData>();
+            std::shared_ptr<CollisionData>      colData   = std::make_shared<CollisionData>();
             std::shared_ptr<CollisionDetection> colDetect = makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), colData);
 
             // Setup the handler
@@ -136,11 +137,11 @@ extern std::shared_ptr<ObjectInteractionPair> makeObjectInteractionPair(
     }
     else if (intType == InteractionType::FemObjToCollidingObj_BoneDrilling)
     {
-        if (std::is_base_of<FeDeformableObject, ObjectType1>::value &&
-            std::is_base_of<CollidingObject, ObjectType2>::value)
+        if (std::is_base_of<FeDeformableObject, ObjectType1>::value
+            && std::is_base_of<CollidingObject, ObjectType2>::value)
         {
             // Setup CD, and collision data
-            std::shared_ptr<CollisionData> colData = std::make_shared<CollisionData>();
+            std::shared_ptr<CollisionData>      colData   = std::make_shared<CollisionData>();
             std::shared_ptr<CollisionDetection> colDetect = makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), colData);
 
             // Setup the handler
diff --git a/Source/Scene/imstkPbdObjectCollisionPair.cpp b/Source/Scene/imstkPbdObjectCollisionPair.cpp
index 14457ec4ba26b6a77618d0019ce0d178426a4f94..9aae06367e46130ad85c211453f9bb8a13c8aa2d 100644
--- a/Source/Scene/imstkPbdObjectCollisionPair.cpp
+++ b/Source/Scene/imstkPbdObjectCollisionPair.cpp
@@ -33,68 +33,65 @@ namespace imstk
 {
 // Pbd Collision will be tested before any step of pbd, then resolved after the solve steps of the two objects
 PbdObjectCollisionPair::PbdObjectCollisionPair(std::shared_ptr<PbdObject> obj1, std::shared_ptr<PbdObject> obj2,
-	CollisionDetection::Type cdType) : CollisionPair(obj1, obj2)
+                                               CollisionDetection::Type cdType) : CollisionPair(obj1, obj2)
 {
-	std::shared_ptr<PbdModel> pbdModel1 = obj1->getPbdModel();
-	std::shared_ptr<PbdModel> pbdModel2 = obj2->getPbdModel();
+    std::shared_ptr<PbdModel> pbdModel1 = obj1->getPbdModel();
+    std::shared_ptr<PbdModel> pbdModel2 = obj2->getPbdModel();
 
+    // Define where collision interaction happens
+    m_computeNodeInputs.first.push_back(pbdModel1->getUpdateCollisionGeometryNode());
+    m_computeNodeInputs.second.push_back(pbdModel2->getUpdateCollisionGeometryNode());
 
-	// Define where collision interaction happens
-	m_computeNodeInputs.first.push_back(pbdModel1->getUpdateCollisionGeometryNode());
-	m_computeNodeInputs.second.push_back(pbdModel2->getUpdateCollisionGeometryNode());
+    m_computeNodeOutputs.first.push_back(pbdModel1->getSolveNode());
+    m_computeNodeOutputs.second.push_back(pbdModel2->getSolveNode());
 
-	m_computeNodeOutputs.first.push_back(pbdModel1->getSolveNode());
-	m_computeNodeOutputs.second.push_back(pbdModel2->getSolveNode());
+    // Define where solver interaction happens
+    m_solveNodeInputs.first.push_back(pbdModel1->getSolveNode());
+    m_solveNodeInputs.second.push_back(pbdModel2->getSolveNode());
 
+    m_solveNodeOutputs.first.push_back(pbdModel1->getUpdateVelocityNode());
+    m_solveNodeOutputs.second.push_back(pbdModel2->getUpdateVelocityNode());
 
-	// Define where solver interaction happens
-	m_solveNodeInputs.first.push_back(pbdModel1->getSolveNode());
-	m_solveNodeInputs.second.push_back(pbdModel2->getSolveNode());
+    // Setup the CD
+    m_colData = std::make_shared<CollisionData>();
+    setCollisionDetection(makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), m_colData));
 
-	m_solveNodeOutputs.first.push_back(pbdModel1->getUpdateVelocityNode());
-	m_solveNodeOutputs.second.push_back(pbdModel2->getUpdateVelocityNode());
+    // Setup the handler
+    std::shared_ptr<PBDCollisionHandling> ch = std::make_shared<PBDCollisionHandling>(CollisionHandling::Side::AB, m_colData, obj1, obj2);
+    setCollisionHandlingAB(ch);
 
-
-	// Setup the CD
-	m_colData = std::make_shared<CollisionData>();
-	setCollisionDetection(makeCollisionDetectionObject(cdType, obj1->getCollidingGeometry(), obj2->getCollidingGeometry(), m_colData));
-
-	// Setup the handler
-	std::shared_ptr<PBDCollisionHandling> ch = std::make_shared<PBDCollisionHandling>(CollisionHandling::Side::AB, m_colData, obj1, obj2);
-	setCollisionHandlingAB(ch);
-
-	// Setup compute node for collision solver (true/critical node)
-	m_collisionSolveNode = std::make_shared<ComputeNode>([ch]() { ch->getCollisionSolver()->solve(); },
-		obj1->getName() + "_vs_" + obj2->getName() + "_CollisionSolver", true);
+    // Setup compute node for collision solver (true/critical node)
+    m_collisionSolveNode = std::make_shared<ComputeNode>([ch]() { ch->getCollisionSolver()->solve(); },
+                obj1->getName() + "_vs_" + obj2->getName() + "_CollisionSolver", true);
 }
 
 void
 PbdObjectCollisionPair::modifyComputeGraph()
 {
-	// Add the collision interaction
-	CollisionPair::modifyComputeGraph();
-
-	// Add a secondary interaction
-	m_objects.first->getComputeGraph()->addNode(m_collisionSolveNode);
-	m_objects.second->getComputeGraph()->addNode(m_collisionSolveNode);
-
-	// Add the solver interaction
-	for (size_t i = 0; i < m_solveNodeInputs.first.size(); i++)
-	{
-		m_objects.first->getComputeGraph()->addEdge(m_solveNodeInputs.first[i], m_collisionSolveNode);
-	}
-	for (size_t i = 0; i < m_solveNodeInputs.second.size(); i++)
-	{
-		m_objects.second->getComputeGraph()->addEdge(m_solveNodeInputs.second[i], m_collisionSolveNode);
-	}
-
-	for (size_t i = 0; i < m_solveNodeOutputs.first.size(); i++)
-	{
-		m_objects.first->getComputeGraph()->addEdge(m_collisionSolveNode, m_solveNodeOutputs.first[i]);
-	}
-	for (size_t i = 0; i < m_solveNodeOutputs.second.size(); i++)
-	{
-		m_objects.second->getComputeGraph()->addEdge(m_collisionSolveNode, m_solveNodeOutputs.second[i]);
-	}
+    // Add the collision interaction
+    CollisionPair::modifyComputeGraph();
+
+    // Add a secondary interaction
+    m_objects.first->getComputeGraph()->addNode(m_collisionSolveNode);
+    m_objects.second->getComputeGraph()->addNode(m_collisionSolveNode);
+
+    // Add the solver interaction
+    for (size_t i = 0; i < m_solveNodeInputs.first.size(); i++)
+    {
+        m_objects.first->getComputeGraph()->addEdge(m_solveNodeInputs.first[i], m_collisionSolveNode);
+    }
+    for (size_t i = 0; i < m_solveNodeInputs.second.size(); i++)
+    {
+        m_objects.second->getComputeGraph()->addEdge(m_solveNodeInputs.second[i], m_collisionSolveNode);
+    }
+
+    for (size_t i = 0; i < m_solveNodeOutputs.first.size(); i++)
+    {
+        m_objects.first->getComputeGraph()->addEdge(m_collisionSolveNode, m_solveNodeOutputs.first[i]);
+    }
+    for (size_t i = 0; i < m_solveNodeOutputs.second.size(); i++)
+    {
+        m_objects.second->getComputeGraph()->addEdge(m_collisionSolveNode, m_solveNodeOutputs.second[i]);
+    }
 }
 }
\ No newline at end of file
diff --git a/Source/Scene/imstkPbdObjectCollisionPair.h b/Source/Scene/imstkPbdObjectCollisionPair.h
index 0a45ac7cae8e9a7b4a9ab9fd81d2018b1281af61..3f19a0ecb659837cab3a998f7504988cc41ca8a6 100644
--- a/Source/Scene/imstkPbdObjectCollisionPair.h
+++ b/Source/Scene/imstkPbdObjectCollisionPair.h
@@ -38,13 +38,13 @@ class PbdObjectCollisionPair : public CollisionPair
 {
 public:
     PbdObjectCollisionPair(std::shared_ptr<PbdObject> obj1, std::shared_ptr<PbdObject> obj2,
-        CollisionDetection::Type cdType = CollisionDetection::Type::MeshToMeshBruteForce);
+                           CollisionDetection::Type cdType = CollisionDetection::Type::MeshToMeshBruteForce);
 
     void modifyComputeGraph() override;
 
 private:
     // Pbd defines two interactions (one at CD and one at solver)
-    Inputs m_solveNodeInputs;
+    Inputs  m_solveNodeInputs;
     Outputs m_solveNodeOutputs;
     std::shared_ptr<ComputeNode> m_collisionSolveNode = nullptr;
 };
diff --git a/Source/Scene/imstkScene.cpp b/Source/Scene/imstkScene.cpp
index 30e8185ad1031a2e0303f0f6e14bded33e1a1e01..f9e3c5fef16ab665d1d243950238c923e96e876d 100644
--- a/Source/Scene/imstkScene.cpp
+++ b/Source/Scene/imstkScene.cpp
@@ -469,11 +469,9 @@ Scene::advance(const double dt)
 
     CollisionDetection::updateInternalOctreeAndDetectCollision();
 
-
     // Execute the computational graph
     m_computeGraphController->execute();
 
-
     // Apply updated forces on device
     for (auto controller : this->getSceneObjectControllers())
     {
diff --git a/Source/Scene/imstkScene.h b/Source/Scene/imstkScene.h
index fed6b0d9abba3ef3189e469299d7e96bb8c685ec..8ccc61507456e99cde4c3b768739c88e4ecb5bb5 100644
--- a/Source/Scene/imstkScene.h
+++ b/Source/Scene/imstkScene.h
@@ -289,12 +289,12 @@ protected:
 
     std::shared_ptr<ComputeGraph> m_computeGraph = nullptr;                      ///> Computational graph
     std::shared_ptr<ComputeGraphController> m_computeGraphController = nullptr;  ///> Controller for the computational graph
-    std::function<void(Scene*)> m_postComputeGraphConfigureCallback = nullptr;
+    std::function<void(Scene*)> m_postComputeGraphConfigureCallback  = nullptr;
 
     std::shared_ptr<ParallelUtils::SpinLock> benchmarkLock = nullptr;
-    std::unordered_map<std::string, double> m_nodeNamesToElapsedTimes; ///> Map of ComputeNode names to elapsed times for benchmarking
+    std::unordered_map<std::string, double>  m_nodeNamesToElapsedTimes; ///> Map of ComputeNode names to elapsed times for benchmarking
 
-    double m_fps         = 0.0;
+    double m_fps = 0.0;
     double m_elapsedTime = 0.0;
 
     bool m_isInitialized = false;
diff --git a/Source/Solvers/imstkDirectLinearSolver.cpp b/Source/Solvers/imstkDirectLinearSolver.cpp
index 77b4018ee1db170783bed698d08968edaaf57dce..1b82a4eb8a063a23f78be90caf956db0cad20986 100644
--- a/Source/Solvers/imstkDirectLinearSolver.cpp
+++ b/Source/Solvers/imstkDirectLinearSolver.cpp
@@ -59,14 +59,20 @@ setSystem(std::shared_ptr<LinearSystem<SparseMatrixd>> newSystem)
 void
 DirectLinearSolver<SparseMatrixd>::solve(const Vectord& rhs, Vectord& x)
 {
-    if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
+    if (!m_linearSystem)
+    {
+        LOG(FATAL) << "Linear system has not been set";
+    }
     x = m_solver.solve(rhs);
 }
 
 void
 DirectLinearSolver<SparseMatrixd>::solve(Vectord& x)
 {
-    if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
+    if (!m_linearSystem)
+    {
+        LOG(FATAL) << "Linear system has not been set";
+    }
     x.setZero();
 
     auto b = m_linearSystem->getRHSVector();
@@ -76,14 +82,20 @@ DirectLinearSolver<SparseMatrixd>::solve(Vectord& x)
 void
 DirectLinearSolver<Matrixd>::solve(const Vectord& rhs, Vectord& x)
 {
-    if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
+    if (!m_linearSystem)
+    {
+        LOG(FATAL) << "Linear system has not been set";
+    }
     x = m_solver.solve(rhs);
 }
 
 void
 DirectLinearSolver<Matrixd>::solve(Vectord& x)
 {
-    if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
+    if (!m_linearSystem)
+    {
+        LOG(FATAL) << "Linear system has not been set";
+    }
     x.setZero();
 
     auto b = m_linearSystem->getRHSVector();
diff --git a/Source/Solvers/imstkDirectLinearSolver.h b/Source/Solvers/imstkDirectLinearSolver.h
index e586d811f03ec1e8fb0b21dd0cdc961d23616017..baaf133e5b79247b15e27777415aaf6e0dd0461e 100644
--- a/Source/Solvers/imstkDirectLinearSolver.h
+++ b/Source/Solvers/imstkDirectLinearSolver.h
@@ -47,7 +47,7 @@ public:
     /// \brief Default constructor/destructor.
     ///
     // DirectLinearSolver() = delete;
-    DirectLinearSolver()  {}
+    DirectLinearSolver() {}
     ~DirectLinearSolver() {};
 
     ///
diff --git a/Source/Solvers/imstkNewtonSolver.cpp b/Source/Solvers/imstkNewtonSolver.cpp
index 3db941c054ff359c6b88d4a1a0da45dec6ae5ebc..28bf71a30b630c900b32351805baa8206d6022eb 100644
--- a/Source/Solvers/imstkNewtonSolver.cpp
+++ b/Source/Solvers/imstkNewtonSolver.cpp
@@ -30,7 +30,7 @@
 
 namespace imstk
 {
-template <> 
+template<>
 NewtonSolver<SparseMatrixd>::NewtonSolver() :
     m_linearSolver(std::make_shared<ConjugateGradient>()),
     m_forcingTerm(0.9),
@@ -43,7 +43,7 @@ NewtonSolver<SparseMatrixd>::NewtonSolver() :
 {
 }
 
-template <>
+template<>
 NewtonSolver<Matrixd>::NewtonSolver() :
     m_linearSolver(std::make_shared<DirectLinearSolver<Matrixd>>()),
     m_forcingTerm(0.9),
@@ -56,8 +56,7 @@ NewtonSolver<Matrixd>::NewtonSolver() :
 {
 }
 
-
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 void
 NewtonSolver<SystemMatrix>::solveGivenState(Vectord& x)
 {
@@ -101,7 +100,7 @@ NewtonSolver<SystemMatrix>::solveGivenState(Vectord& x)
     }
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 void
 NewtonSolver<SystemMatrix>::solve()
 {
@@ -112,14 +111,14 @@ NewtonSolver<SystemMatrix>::solve()
     }
 
     size_t      iterNum;
-    const auto& u  = this->m_nonLinearSystem->getUnknownVector();
-    Vectord     du = u; // make this a class member in future
-    double error0 = MAX_D;
+    const auto& u      = this->m_nonLinearSystem->getUnknownVector();
+    Vectord     du     = u; // make this a class member in future
+    double      error0 = MAX_D;
 
     double epsilon = m_relativeTolerance * m_relativeTolerance;
     for (iterNum = 0; iterNum < m_maxIterations; ++iterNum)
     {
-        double error  = updateJacobian(u);
+        double error = updateJacobian(u);
 
         if (iterNum == 0)
         {
@@ -144,7 +143,7 @@ NewtonSolver<SystemMatrix>::solve()
     }
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 double
 NewtonSolver<SystemMatrix>::updateJacobian(const Vectord& x)
 {
@@ -171,7 +170,7 @@ NewtonSolver<SystemMatrix>::updateJacobian(const Vectord& x)
     return std::sqrt(b.dot(b));
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 void
 NewtonSolver<SystemMatrix>::updateForcingTerm(const double ratio, const double stopTolerance, const double fnorm)
 {
@@ -188,29 +187,28 @@ NewtonSolver<SystemMatrix>::updateForcingTerm(const double ratio, const double s
     m_forcingTerm = std::max(std::min(eta, m_etaMax), 0.5 * stopTolerance / fnorm);
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 // std::shared_ptr<NewtonSolver<SystemMatrix>::LinearSolverType>
 auto
-NewtonSolver<SystemMatrix>::getLinearSolver() const -> std::shared_ptr<LinearSolverType>
+NewtonSolver<SystemMatrix>::getLinearSolver() const->std::shared_ptr<LinearSolverType>
 {
     return m_linearSolver;
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 void
 NewtonSolver<SystemMatrix>::setAbsoluteTolerance(const double aTolerance)
 {
     m_absoluteTolerance = aTolerance;
 }
 
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 double
 NewtonSolver<SystemMatrix>::getAbsoluteTolerance() const
 {
     return m_absoluteTolerance;
 }
 
-
 template class NewtonSolver<SparseMatrixd>;
 template class NewtonSolver<Matrixd>;
 } // imstk
diff --git a/Source/Solvers/imstkNewtonSolver.h b/Source/Solvers/imstkNewtonSolver.h
index 47af61a2a0d09233b9eb99cb30ba59ce0e1c2617..2dded1606af66aec5fecab42a955738777bba738 100644
--- a/Source/Solvers/imstkNewtonSolver.h
+++ b/Source/Solvers/imstkNewtonSolver.h
@@ -33,7 +33,7 @@ namespace imstk
 ///     root of F or fails. Global convergence is achieved using a line search sub-process
 ///     and the Armijo rule
 ///
-template <typename SystemMatrix> 
+template<typename SystemMatrix>
 class NewtonSolver : public NonLinearSolver<SystemMatrix>
 {
 public:
@@ -75,9 +75,9 @@ public:
     /// \param newLinearSolver Linear solver pointer
     ///
     void setLinearSolver(std::shared_ptr<LinearSolverType> newLinearSolver)
-	{
-		m_linearSolver = newLinearSolver;
-	}
+    {
+        m_linearSolver = newLinearSolver;
+    }
 
     ///
     /// \brief Get LinearSolver
@@ -207,7 +207,7 @@ public:
     virtual void setToSemiImplicit() override
     {
         this->m_isSemiImplicit = true;
-        m_maxIterations  = 1;
+        m_maxIterations = 1;
     }
 
 private:
diff --git a/Source/Solvers/imstkNonLinearSolver.cpp b/Source/Solvers/imstkNonLinearSolver.cpp
index bf7dcc7cbe2675d4ad8ae48fffde8412749d76c3..6eac4676c94c171ba1a6e480f8ee1160401c8ffa 100644
--- a/Source/Solvers/imstkNonLinearSolver.cpp
+++ b/Source/Solvers/imstkNonLinearSolver.cpp
@@ -24,7 +24,7 @@
 
 namespace imstk
 {
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 NonLinearSolver<SystemMatrix>::NonLinearSolver() : m_sigma(std::array<double, 2>
     {
         { 0.1, 0.5 }
@@ -38,7 +38,7 @@ NonLinearSolver<SystemMatrix>::NonLinearSolver() : m_sigma(std::array<double, 2>
                       };
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 double
 NonLinearSolver<SystemMatrix>::armijo(const Vectord& dx, Vectord& x, const double previousFnorm)
 {
@@ -99,7 +99,7 @@ NonLinearSolver<SystemMatrix>::armijo(const Vectord& dx, Vectord& x, const doubl
     return currentFnorm;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 void
 NonLinearSolver<SystemMatrix>::parabolicModel(const std::array<double, 3>& fnorm, std::array<double, 3>& lambda)
 {
@@ -134,56 +134,56 @@ NonLinearSolver<SystemMatrix>::parabolicModel(const std::array<double, 3>& fnorm
     lambda[0] = newLambda;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 void
 NonLinearSolver<SystemMatrix>::setSigma(const std::array<double, 2>& newSigma)
 {
     m_sigma = newSigma;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 const std::array<double, 2>&
 NonLinearSolver<SystemMatrix>::getSigma() const
 {
     return m_sigma;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 void
 NonLinearSolver<SystemMatrix>::setAlpha(const double newAlpha)
 {
     m_alpha = newAlpha;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 double
 NonLinearSolver<SystemMatrix>::getAlpha() const
 {
     return m_alpha;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 void
 NonLinearSolver<SystemMatrix>::setArmijoMax(const size_t newArmijoMax)
 {
     m_armijoMax = newArmijoMax;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 size_t
 NonLinearSolver<SystemMatrix>::getArmijoMax() const
 {
     return m_armijoMax;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 void
 NonLinearSolver<SystemMatrix>::setSystem(std::shared_ptr<NonLinearSystem<SystemMatrix>> newSystem)
 {
     m_nonLinearSystem = newSystem;
 }
 
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 std::shared_ptr<NonLinearSystem<SystemMatrix>>
 NonLinearSolver<SystemMatrix>::getSystem() const
 {
diff --git a/Source/Solvers/imstkNonLinearSolver.h b/Source/Solvers/imstkNonLinearSolver.h
index 982ecb2c794471d0ed3a9a2c92c89965325b0406..f576ea458b2376f5c8574c01976adee1bf2af6ef 100644
--- a/Source/Solvers/imstkNonLinearSolver.h
+++ b/Source/Solvers/imstkNonLinearSolver.h
@@ -35,7 +35,7 @@ namespace imstk
 ///
 /// \brief Base class for non-linear solvers
 ///
-template <typename SystemMatrix>
+template<typename SystemMatrix>
 class NonLinearSolver : public SolverBase
 {
 public:
@@ -115,9 +115,9 @@ public:
     /// \param newUpdateIterate Function used to update iterates. Default: x+=dx.
     ///
     void setUpdateIterate(const UpdateIterateType& newUpdateIterate)
-	{
-		m_updateIterate = newUpdateIterate;
-	}
+    {
+        m_updateIterate = newUpdateIterate;
+    }
 
     ///
     /// \brief Set the Newton solver to be fully implicit
@@ -136,12 +136,12 @@ public:
     }
 
 protected:
-    std::array<double, 2> m_sigma;                      ///< Safeguarding bounds for the line search
-    double m_alpha;                                     ///< Parameter to measure decrease
-    size_t m_armijoMax;                                 ///< Maximum number of step length reductions
+    std::array<double, 2> m_sigma;                                    ///< Safeguarding bounds for the line search
+    double m_alpha;                                                   ///< Parameter to measure decrease
+    size_t m_armijoMax;                                               ///< Maximum number of step length reductions
 
     std::shared_ptr<NonLinearSystem<SystemMatrix>> m_nonLinearSystem; ///< System of non-linear equations
-    UpdateIterateType m_updateIterate;                  ///< Update iteration function
-    bool m_isSemiImplicit = true;                       ///> Semi-Implicit solver
+    UpdateIterateType m_updateIterate;                                ///< Update iteration function
+    bool m_isSemiImplicit = true;                                     ///> Semi-Implicit solver
 };
 } // imstk
diff --git a/Source/Solvers/imstkNonLinearSystem.cpp b/Source/Solvers/imstkNonLinearSystem.cpp
index 92c688c54b8a1ac4b4e0409814c7fdfb58b6bc80..0ce7bc971b8f0d85199110ca511ba94aad6d975f 100644
--- a/Source/Solvers/imstkNonLinearSystem.cpp
+++ b/Source/Solvers/imstkNonLinearSystem.cpp
@@ -26,33 +26,33 @@
 
 namespace imstk
 {
-template <typename Matrix>
+template<typename Matrix>
 NonLinearSystem<Matrix>::NonLinearSystem(const VectorFunctionType& F, const MatrixFunctionType& dF) : m_F(F), m_dF(dF)
 {
 }
 
-template <typename Matrix>
+template<typename Matrix>
 void
 NonLinearSystem<Matrix>::setFunction(const VectorFunctionType& function)
 {
     m_F = function;
 }
 
-template <typename Matrix>
+template<typename Matrix>
 void
 NonLinearSystem<Matrix>::setJacobian(const MatrixFunctionType& function)
 {
     m_dF = function;
 }
 
-template <typename Matrix>
+template<typename Matrix>
 const Vectord&
 NonLinearSystem<Matrix>::evaluateF(const Vectord& x, const bool isSemiImplicit)
 {
     return m_F(x, isSemiImplicit);
 }
 
-template <typename Matrix>
+template<typename Matrix>
 const Matrix&
 NonLinearSystem<Matrix>::evaluateJacobian(const Vectord& x)
 {
@@ -61,5 +61,4 @@ NonLinearSystem<Matrix>::evaluateJacobian(const Vectord& x)
 
 template class NonLinearSystem<SparseMatrixd>;
 template class NonLinearSystem<Matrixd>;
-
 } //imstk
diff --git a/Source/Solvers/imstkNonLinearSystem.h b/Source/Solvers/imstkNonLinearSystem.h
index d631063b9b287bb562a671059b09ea35e92ca1bf..2a49cc24362b31b3c083113b11f83c76d012bac7 100644
--- a/Source/Solvers/imstkNonLinearSystem.h
+++ b/Source/Solvers/imstkNonLinearSystem.h
@@ -31,7 +31,7 @@ namespace imstk
 ///
 /// \brief Base class for a multi-variable nonlinear system
 ///
-template <typename Matrix>
+template<typename Matrix>
 class NonLinearSystem
 {
 public:
diff --git a/Source/Solvers/imstkPbdSolver.cpp b/Source/Solvers/imstkPbdSolver.cpp
index 16c5eb09e2b8a89a3980aaf54281bc525a23aa3d..25accfa1988a429ad1ad5d9188014056f7535dd6 100644
--- a/Source/Solvers/imstkPbdSolver.cpp
+++ b/Source/Solvers/imstkPbdSolver.cpp
@@ -69,7 +69,6 @@ PbdSolver::solve()
     }
 }
 
-
 PbdCollisionSolver::PbdCollisionSolver() :
     m_collisionConstraints(std::make_shared<std::list<PBDCollisionConstraintVector*>>()),
     m_collisionConstraintsData(std::make_shared<std::list<CollisionConstraintData>>())
@@ -78,8 +77,8 @@ PbdCollisionSolver::PbdCollisionSolver() :
 
 void
 PbdCollisionSolver::addCollisionConstraints(PBDCollisionConstraintVector* constraints,
-    std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
-    std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB)
+                                            std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
+                                            std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB)
 {
     m_collisionConstraints->push_back(constraints);
     m_collisionConstraintsData->push_back({ posA, invMassA, posB, invMassB });
@@ -97,11 +96,11 @@ PbdCollisionSolver::solve()
             std::list<CollisionConstraintData>::iterator colDataIter = m_collisionConstraintsData->begin();
             for (auto constraintList : *m_collisionConstraints)
             {
-                CollisionConstraintData colData  = *colDataIter;
-                StdVectorOfVec3d&       posA     = *colData.m_posA;
-                const StdVectorOfReal&  invMassA = *colData.m_invMassA;
-                StdVectorOfVec3d&       posB     = *colData.m_posB;
-                const StdVectorOfReal&  invMassB = *colData.m_invMassB;
+                CollisionConstraintData             colData     = *colDataIter;
+                StdVectorOfVec3d&                   posA        = *colData.m_posA;
+                const StdVectorOfReal&              invMassA    = *colData.m_invMassA;
+                StdVectorOfVec3d&                   posB        = *colData.m_posB;
+                const StdVectorOfReal&              invMassB    = *colData.m_invMassB;
                 const PBDCollisionConstraintVector& constraints = *constraintList;
                 for (size_t j = 0; j < constraints.size(); j++)
                 {
@@ -114,5 +113,4 @@ PbdCollisionSolver::solve()
         m_collisionConstraintsData->clear();
     }
 }
-
 } // end namespace imstk
diff --git a/Source/Solvers/imstkPbdSolver.h b/Source/Solvers/imstkPbdSolver.h
index 87c3cd9124286aec3c691d3a0e4e42cea77222f8..0622a90ec0911e88eed660e2da4c17ddce234826 100644
--- a/Source/Solvers/imstkPbdSolver.h
+++ b/Source/Solvers/imstkPbdSolver.h
@@ -30,9 +30,9 @@ namespace imstk
 struct CollisionConstraintData
 {
     CollisionConstraintData(std::shared_ptr<StdVectorOfVec3d> posA,
-        std::shared_ptr<StdVectorOfReal> invMassA,
-        std::shared_ptr<StdVectorOfVec3d> posB,
-        std::shared_ptr<StdVectorOfReal> invMassB) :
+                            std::shared_ptr<StdVectorOfReal>  invMassA,
+                            std::shared_ptr<StdVectorOfVec3d> posB,
+                            std::shared_ptr<StdVectorOfReal>  invMassB) :
         m_posA(posA), m_invMassA(invMassA), m_posB(posB), m_invMassB(invMassB)
     {
     }
@@ -98,7 +98,7 @@ public:
     void solve() override;
 
 private:
-    size_t m_iterations = 20; ///> Number of NL Gauss-Seidel iterations for regular constraints
+    size_t m_iterations = 20;                                                             ///> Number of NL Gauss-Seidel iterations for regular constraints
 
     std::shared_ptr<std::vector<PBDConstraintVector>> m_partitionedConstraints = nullptr; ///> Set of vector'd/partitioned pbd constraints
     std::shared_ptr<PBDConstraintVector> m_constraints = nullptr;                         ///> Vector of constraints
@@ -123,8 +123,8 @@ public:
     /// \brief Add the global collision contraints to this solver
     ///
     void addCollisionConstraints(PBDCollisionConstraintVector* constraints,
-        std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
-        std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB);
+                                 std::shared_ptr<StdVectorOfVec3d> posA, std::shared_ptr<StdVectorOfReal> invMassA,
+                                 std::shared_ptr<StdVectorOfVec3d> posB, std::shared_ptr<StdVectorOfReal> invMassB);
 
     ///
     /// \brief Solve the non linear system of equations G(x)=0 using Newton's method.
@@ -132,7 +132,7 @@ public:
     void solve() override;
 
 private:
-    size_t m_collisionIterations = 5; ///> Number of NL Gauss-Seidel iterations for collision constraints
+    size_t m_collisionIterations = 5;                                                               ///> Number of NL Gauss-Seidel iterations for collision constraints
 
     std::shared_ptr<std::list<PBDCollisionConstraintVector*>> m_collisionConstraints     = nullptr; ///< Collision contraints charged to this solver
     std::shared_ptr<std::list<CollisionConstraintData>>       m_collisionConstraintsData = nullptr;
diff --git a/Source/apiUtilities/imstkAPIUtilities.h b/Source/apiUtilities/imstkAPIUtilities.h
index f39515b79aa620b9c7a23580a178117a2fe608cc..8a929f125554ba9c91a75f12b5e0320696cfa25e 100644
--- a/Source/apiUtilities/imstkAPIUtilities.h
+++ b/Source/apiUtilities/imstkAPIUtilities.h
@@ -33,7 +33,7 @@ class PointSet;
 class SceneManager;
 class Scene;
 class FEMDeformableBodyModel;
-template <typename Matrix>
+template<typename Matrix>
 class NonLinearSystem;
 class CollidingObject;
 class Graph;