diff --git a/Collision/PlaneToMeshCollision.cpp b/Collision/PlaneToMeshCollision.cpp
index 0b13d3f049174296af10a127962ef5d338d2b993..2799d8770dcaff24074dca13cd057b08866eed12 100644
--- a/Collision/PlaneToMeshCollision.cpp
+++ b/Collision/PlaneToMeshCollision.cpp
@@ -34,30 +34,28 @@
 
 void PlaneToMeshCollision::doComputeCollision(std::shared_ptr<CollisionPair> pair)
 {
-    auto mesh = std::static_pointer_cast<MeshCollisionModel>(pair->getFirst());
-    auto plane = std::static_pointer_cast<PlaneCollisionModel>(pair->getSecond());
+    auto meshModel = std::static_pointer_cast<MeshCollisionModel>(pair->getFirst());
+    auto planeModel = std::static_pointer_cast<PlaneCollisionModel>(pair->getSecond());
 
-    if (!mesh || !plane)
+    if (!meshModel || !planeModel)
     {
         return;
     }
 
-    double d;
-    core::Vec3d planeNormal = plane->getPlaneModel()->getUnitNormal();
-    core::Vec3d planePos = plane->getPlaneModel()->getPoint();
+    core::Vec3d normal = planeModel->getPlaneModel()->getUnitNormal();
+    core::Vec3d contactPoint = planeModel->getPlaneModel()->getPoint();
 
-    core::Vec3d vert;
     pair->clearContacts();
-    for (size_t i = 0; i < mesh->getVertices().size(); i++)//const auto& vertex : mesh->getVertices()
+    const auto &vertices = meshModel->getVertices();
+    for (size_t i = 0, end = vertices.size(); i < end; ++i)
     {
-        vert = mesh->getVertices()[i];
-        d = planeNormal.dot(vert - planePos);
+        auto vertex = vertices[i];
+        auto d = normal.dot(vertex - contactPoint);
 
         if (d < std::numeric_limits<float>::epsilon())
         {
-            pair->addContact(d, vert, i, planeNormal);// Create contact
+            pair->addContact(meshModel, d, vertex, i, normal);// Create contact
+            pair->addContact(planeModel, d, contactPoint, i, -normal);// Create contact
         }
     }
-    /*std::cout << "@ Collision detection\n";
-    pair->printCollisionPairs();*/
 }
diff --git a/ContactHandling/PenaltyContactFemToStatic.cpp b/ContactHandling/PenaltyContactFemToStatic.cpp
index 04b09d775fc98d2f5fa2279dd532072ffb0fdc7d..4fab3cdbc78700f06a0c73ed367c528f5ce71b88 100644
--- a/ContactHandling/PenaltyContactFemToStatic.cpp
+++ b/ContactHandling/PenaltyContactFemToStatic.cpp
@@ -44,35 +44,8 @@ PenaltyContactFemToStatic::~PenaltyContactFemToStatic()
 
 void PenaltyContactFemToStatic::computeUnilateralContactForces()
 {
-    int nodeDofID;
-    const double stiffness = 1.0e4, damping = 1.0e5;
-    core::Vec3d velocityProjection;
-
-    std::vector<std::shared_ptr<Contact>> contactInfo = this->getCollisionPairs()->getContacts();
-
-    if (this->getSecondSceneObject()->getType() == core::ClassType::VegaFemSceneObject
-        && this->getFirstSceneObject()->getType() == core::ClassType::StaticSceneObject)
-    {
-        auto femSceneObject = std::static_pointer_cast<VegaFemSceneObject>(this->getSecondSceneObject());
-
-        femSceneObject->setContactForcesToZero();
-        core::Vec3d force;
-        for (size_t i = 0; i < contactInfo.size(); i++)
-        {
-            nodeDofID = 3 * contactInfo[i]->index;
-            velocityProjection = femSceneObject->getVelocityOfNodeWithDofID(nodeDofID);
-            velocityProjection = contactInfo[i]->normal.dot(velocityProjection) * contactInfo[i]->normal;
-
-            force = stiffness * -contactInfo[i]->depth * contactInfo[i]->normal - damping * velocityProjection;
-
-            femSceneObject->setContactForceOfNodeWithDofID(nodeDofID, contactInfo[i]->point, force);
-
-        }
-    }
-    else
-    {
-        std::cout << "Error: Scene objects don't match the required types in 'PenaltyContactFemToStatic' \n";
-    }
+    auto femSceneObject = std::static_pointer_cast<VegaFemSceneObject>(this->getSecondSceneObject());
+    this->computeForces(femSceneObject);
 }
 
 void PenaltyContactFemToStatic::computeBilateralContactForces()
diff --git a/ContactHandling/PenaltyContactHandling.cpp b/ContactHandling/PenaltyContactHandling.cpp
index d6e0e291d3dd052a74f4e5c5066300f5e6e91714..3caa1d494550b4c4d22684dc3e7844bb76d770ab 100644
--- a/ContactHandling/PenaltyContactHandling.cpp
+++ b/ContactHandling/PenaltyContactHandling.cpp
@@ -25,14 +25,20 @@
 #include "Simulators/VegaFemSceneObject.h"
 #include "Core/CollisionPair.h"
 
-PenaltyContactHandling::PenaltyContactHandling(bool typeBilateral) : ContactHandling(typeBilateral)
+PenaltyContactHandling::PenaltyContactHandling(bool typeBilateral) :
+    ContactHandling(typeBilateral),
+    stiffness(1e4),
+    damping(1e5)
+
 {
 }
 
-PenaltyContactHandling::PenaltyContactHandling( bool typeBilateral,
-                                                    const std::shared_ptr<SceneObject>& sceneObjFirst,
-                                                    const std::shared_ptr<SceneObject>& sceneObjSecond)
-                                                    : ContactHandling(typeBilateral,sceneObjFirst,sceneObjSecond)
+PenaltyContactHandling::PenaltyContactHandling(bool typeBilateral,
+                                               const std::shared_ptr<SceneObject>& sceneObjFirst,
+                                               const std::shared_ptr<SceneObject>& sceneObjSecond) :
+    ContactHandling(typeBilateral,sceneObjFirst,sceneObjSecond),
+    stiffness(1e4),
+    damping(1e5)
 {
 }
 
diff --git a/ContactHandling/PenaltyContactHandling.h b/ContactHandling/PenaltyContactHandling.h
index 1ac597166d1dbb51a4a3c52288794525e5cb2407..444f61ebd064a4a4f9a12b861b82b6e5445f1271 100644
--- a/ContactHandling/PenaltyContactHandling.h
+++ b/ContactHandling/PenaltyContactHandling.h
@@ -26,6 +26,8 @@
 // SimMedTK includes
 #include "Core/ContactHandling.h"
 
+#include "Core/CollisionPair.h"
+
 /// \brief Penalty based for contact handling
 class PenaltyContactHandling : public ContactHandling
 {
@@ -33,8 +35,8 @@ public:
     PenaltyContactHandling(bool typeBilateral);
 
     PenaltyContactHandling(bool typeBilateral,
-                             const std::shared_ptr<SceneObject>& sceneObjFirst,
-                             const std::shared_ptr<SceneObject>& sceneObjSecond);
+                           const std::shared_ptr<SceneObject>& sceneObjFirst,
+                           const std::shared_ptr<SceneObject>& sceneObjSecond);
 
     virtual ~PenaltyContactHandling();
 
@@ -45,6 +47,73 @@ public:
 
     /// \brief Get the forces on both the scene objects using penalty method
     virtual void computeBilateralContactForces() = 0;
+
+    /// \brief Get the forces on both the scene objects using penalty method
+    virtual void computeForces(std::shared_ptr<SceneObject> sceneObject)
+    {
+        if(sceneObject->computeContactForce())
+        {
+            auto model = sceneObject->getModel();
+            if(!model)
+            {
+                return;
+            }
+
+            auto contactInfo = this->getCollisionPairs()->getContacts(model);
+            sceneObject->setContactForcesToZero();
+            core::Vec3d force;
+            core::Vec3d velocityProjection;
+            int nodeDofID;
+
+            for(size_t i = 0; i < contactInfo.size(); ++i)
+            {
+                nodeDofID = 3 * contactInfo[i]->index;
+                velocityProjection = sceneObject->getVelocity(nodeDofID);
+                velocityProjection = contactInfo[i]->normal.dot(velocityProjection) * contactInfo[i]->normal;
+
+                force = -stiffness * contactInfo[i]->depth * contactInfo[i]->normal - damping * velocityProjection;
+
+                sceneObject->setContactForce(nodeDofID, contactInfo[i]->point, force);
+            }
+            return;
+        }
+    }
+
+    ///
+    /// \brief Set stiffness coefficient
+    ///
+    void setStiffness(const double stiffnessCoeff)
+    {
+        this->stiffness = stiffnessCoeff;
+    }
+
+    ///
+    /// \brief Get the stiffness coefficient
+    ///
+    double getStiffness() const
+    {
+        return this->stiffness;
+    }
+
+    ///
+    /// \brief Set stiffness coefficient
+    ///
+    void setDamping(const double dampingValue)
+    {
+        this->damping = dampingValue;
+    }
+
+    ///
+    /// \brief Get the stiffness coefficient
+    ///
+    double getDamping() const
+    {
+        return this->damping;
+    }
+
+private:
+    double stiffness;
+    double damping;
 };
 
 #endif // SMPENALTY_CONTACTHANDLING_H
diff --git a/ContactHandling/UnitTests/PenaltyContactHandlingSpec.cpp b/ContactHandling/UnitTests/PenaltyContactHandlingSpec.cpp
index 42e2e409b9125010dd5006d4173f440f1fa928bd..140cf9d7e763bffd27e58e54a97c6665bca13a10 100644
--- a/ContactHandling/UnitTests/PenaltyContactHandlingSpec.cpp
+++ b/ContactHandling/UnitTests/PenaltyContactHandlingSpec.cpp
@@ -107,7 +107,7 @@ go_bandit([]() {
             f.push_back(0);
             f.push_back(0);
 
-            fem->setContactForceOfNodeWithDofID(0,core::Vec3d(0,0,0));
+            fem->setContactForce(0,core::Vec3d(0,0,0));
             auto &contactForce = fem->getContactForces();
 
             handler->resolveContacts();
diff --git a/Core/CollisionPair.cpp b/Core/CollisionPair.cpp
index f5a7a859e9ab780dc369726f834c936360449fcb..a93d8bccce6227a61f3ac0cfa852f37e2d1681c9 100644
--- a/Core/CollisionPair.cpp
+++ b/Core/CollisionPair.cpp
@@ -22,65 +22,114 @@
 //---------------------------------------------------------------------------
 
 #include "Core/CollisionPair.h"
+#include "Core/Model.h"
 
 CollisionPair::CollisionPair() {}
+
+//---------------------------------------------------------------------------
 CollisionPair::~CollisionPair() {}
 
-void CollisionPair::setModels(const std::shared_ptr< Model >& first,
-                                const std::shared_ptr< Model >& second)
+//---------------------------------------------------------------------------
+void CollisionPair::setModels(std::shared_ptr<Model> first,
+                              std::shared_ptr<Model> second)
 {
     this->modelRepresentations.first = first;
     this->modelRepresentations.second = second;
 }
 
-void CollisionPair::addContact(const double& penetrationDepth, const core::Vec3d& vert, const int index, const core::Vec3d& contactNornmal)
+//---------------------------------------------------------------------------
+void CollisionPair::addContact(const double& penetrationDepth,
+                               const core::Vec3d& vert,
+                               const int index,
+                               const core::Vec3d& contactNornmal)
 {
-    auto contact = std::make_shared<Contact>(penetrationDepth, vert, index, contactNornmal);
-    //std::shared_ptr<Contact> contact(new Contact(penetrationDepth, vert, index, contactNornmal));
+    auto contact = std::make_shared<Contact>(penetrationDepth,
+                                             vert,
+                                             index,
+                                             contactNornmal);
     this->contacts.emplace_back(contact);
 }
 
-const std::pair< std::shared_ptr< Model >, std::shared_ptr< Model > >&
+//---------------------------------------------------------------------------
+void CollisionPair::addContact(std::shared_ptr<Model> model,
+                               const double& penetrationDepth,
+                               const core::Vec3d& vert,
+                               const int index,
+                               const core::Vec3d& contactNornmal)
+{
+    auto contact = std::make_shared<Contact>(penetrationDepth,
+                                             vert,
+                                             index,
+                                             contactNornmal);
+    this->contacts.emplace_back(contact);
+    this->modelContacts[model].emplace_back(contact);
+}
+
+//---------------------------------------------------------------------------
+const std::pair< std::shared_ptr<Model>, std::shared_ptr<Model> >&
 CollisionPair::getModels() const
 {
-    return modelRepresentations;
+    return this->modelRepresentations;
 }
 
+//---------------------------------------------------------------------------
 void CollisionPair::clearContacts()
 {
-    contacts.clear();
+    this->contacts.clear();
 }
 
+//---------------------------------------------------------------------------
 int CollisionPair::getNumberOfContacts()
 {
-    return contacts.size();
+    return this->contacts.size();
 }
 
-std::shared_ptr< Model > CollisionPair::getFirst()
+//---------------------------------------------------------------------------
+std::shared_ptr<Model> CollisionPair::getFirst()
 {
     return this->modelRepresentations.first;
 }
 
-std::shared_ptr< Model > CollisionPair::getSecond()
+//---------------------------------------------------------------------------
+std::shared_ptr<Model> CollisionPair::getSecond()
 {
     return this->modelRepresentations.second;
 }
 
+//---------------------------------------------------------------------------
 bool CollisionPair::hasContacts()
 {
     return !this->contacts.empty();
 }
 
+//---------------------------------------------------------------------------
 std::vector< std::shared_ptr< Contact > >& CollisionPair::getContacts()
 {
-    return contacts;
+    return this->contacts;
 }
 
-const std::vector< std::shared_ptr< Contact > >& CollisionPair::getContacts() const
+//---------------------------------------------------------------------------
+const std::vector< std::shared_ptr< Contact > >&
+CollisionPair::getContacts() const
 {
-    return contacts;
+    return this->contacts;
 }
 
+//---------------------------------------------------------------------------
+std::vector< std::shared_ptr< Contact > >&
+CollisionPair::getContacts(const std::shared_ptr<Model> &model)
+{
+    return this->modelContacts.at(model);
+}
+
+//---------------------------------------------------------------------------
+const std::vector< std::shared_ptr< Contact > >&
+CollisionPair::getContacts(const std::shared_ptr<Model> &model) const
+{
+    return this->contacts;
+}
+
+//---------------------------------------------------------------------------
 void CollisionPair::printCollisionPairs()
 {
     std::cout << "# Contacts: " << this->contacts.size() << std::endl;
diff --git a/Core/CollisionPair.h b/Core/CollisionPair.h
index e50f01de8f19d45fe56ecaba7fde44485c760771..a257fad8af7bb284ccc1442ebb9d371163066fd3 100644
--- a/Core/CollisionPair.h
+++ b/Core/CollisionPair.h
@@ -32,6 +32,7 @@
 #include <memory>
 #include <vector>
 #include <iostream>
+#include <map>
 
 class Model;
 
@@ -45,7 +46,21 @@ public:
              const core::Vec3d& p,
              const int ind,
              const core::Vec3d& contactNornmal) :
-                depth(penetrationDepth), point(p), normal(contactNornmal), index(ind){}
+                depth(penetrationDepth),
+                point(p),
+                normal(contactNornmal),
+                index(ind),
+                model(nullptr){}
+    Contact (std::shared_ptr<Model> m,
+             const double penetrationDepth,
+             const core::Vec3d& p,
+             const int ind,
+             const core::Vec3d& contactNornmal) :
+                depth(penetrationDepth),
+                point(p),
+                normal(contactNornmal),
+                index(ind),
+                model(m){}
 
     void printInfo()
     {
@@ -59,6 +74,7 @@ public:
     core::Vec3d point;
     core::Vec3d normal;
     int index;
+    std::shared_ptr<Model> model;
 };
 
 ///
@@ -74,8 +90,8 @@ public:
     ///
     /// @brief Set the pair of collision models
     ///
-    void setModels(const std::shared_ptr<Model>& first,
-                   const std::shared_ptr<Model>& second );
+    void setModels(std::shared_ptr<Model> first,
+                   std::shared_ptr<Model> second );
 
     ///
     /// @brief Get the pair of collision models
@@ -91,6 +107,15 @@ public:
                      const int index,
                      const core::Vec3d& contactNornmal);
 
+    ///
+    /// @brief Add contact between the models
+    ///
+    void addContact( std::shared_ptr<Model> model,
+                     const double& penetrationDepth,
+                     const core::Vec3d& vert,
+                     const int index,
+                     const core::Vec3d& contactNornmal);
+
     ///
     /// @brief Clear contact list
     ///
@@ -122,6 +147,12 @@ public:
     std::vector<std::shared_ptr<Contact>> &getContacts();
     const std::vector<std::shared_ptr<Contact>> &getContacts() const;
 
+    ///
+    /// @brief Returns contact array for a particular model
+    ///
+    std::vector<std::shared_ptr<Contact>> &getContacts(const std::shared_ptr<Model> &model);
+    const std::vector<std::shared_ptr<Contact>> &getContacts(const std::shared_ptr<Model> &model) const;
+
     ///
     /// @brief Returns contact array for these two models
     ///
@@ -129,9 +160,10 @@ public:
 
 private:
     std::pair<std::shared_ptr<Model>,
-        std::shared_ptr<Model>> modelRepresentations; // Models
-
-    std::vector<std::shared_ptr<Contact>> contacts; // List of contacts
+        std::shared_ptr<Model>> modelRepresentations; //!< Models
+    std::vector<std::shared_ptr<Contact>> contacts; //!< List of contacts
+    std::map<std::shared_ptr<Model>,
+               std::vector<std::shared_ptr<Contact>>> modelContacts; //!< Contacts per model
 };
 
 #endif // SMCOLLISIONPAIR_H
diff --git a/Core/ContactHandling.cpp b/Core/ContactHandling.cpp
index 83f259e13df58c5b0884ce5a605b7700cff41765..bb14610d18b9810a067d5e94411d4c76b5271fba 100644
--- a/Core/ContactHandling.cpp
+++ b/Core/ContactHandling.cpp
@@ -23,6 +23,8 @@
 
 #include "Core/ContactHandling.h"
 
+#include "Core/CollisionPair.h"
+
 ContactHandling::ContactHandling(const bool typeBilateral)
 {
     if (typeBilateral)
@@ -69,12 +71,12 @@ void ContactHandling::setSceneObjects(const std::shared_ptr< SceneObject > first
 
 void ContactHandling::setCollisionPairs(const std::shared_ptr< CollisionPair > colPair)
 {
-    collisionPairs = colPair;
+    collisionPair = colPair;
 }
 
 std::shared_ptr<CollisionPair> ContactHandling::getCollisionPairs() const
 {
-    return collisionPairs;
+    return collisionPair;
 }
 
 ContactHandlingType ContactHandling::getContactHandlingType() const
diff --git a/Core/ContactHandling.h b/Core/ContactHandling.h
index 848f9cb4aa52770008f36e8c2d69712b4fd2247b..31cce97d7389f4801ceb2c491e9dcd705279d9dc 100644
--- a/Core/ContactHandling.h
+++ b/Core/ContactHandling.h
@@ -81,6 +81,7 @@ public:
     /// \brief Implementation of how the contacts between colliding
     /// objects is resolved
     virtual void resolveContacts() = 0;
+
 protected:
 
     ContactHandlingType type;
@@ -89,7 +90,7 @@ protected:
 
     std::pair<std::shared_ptr<SceneObject>, std::shared_ptr<SceneObject>> collidingSceneObjects;
 
-    std::shared_ptr<CollisionPair> collisionPairs;
+    std::shared_ptr<CollisionPair> collisionPair;
 };
 
 #endif //SMCONTACTHANDLING_H
diff --git a/Core/SceneObject.cpp b/Core/SceneObject.cpp
index 812466d562ce65121bc07fdc3f2bbfc73273df54..5aac7010a367eaf7bf529e3e4668eeed747db067 100644
--- a/Core/SceneObject.cpp
+++ b/Core/SceneObject.cpp
@@ -35,6 +35,7 @@ SceneObject::SceneObject()
     flags.isViewerInit = false;
     flags.isSimulatorInit = false;
     name = "SceneObject" + std::to_string(this->getUniqueId()->getId());
+    hasContactForces = false;
 }
 
 SceneObject::~SceneObject()
diff --git a/Core/SceneObject.h b/Core/SceneObject.h
index 9e3b4aa21929f1671157f55ed3f3a17b7adf9e74..e97f95b4b103c2e9a8e32a8401ace83fa0e9bb36 100644
--- a/Core/SceneObject.h
+++ b/Core/SceneObject.h
@@ -24,9 +24,10 @@
 #ifndef SMSCENEOBJECT_H
 #define SMSCENEOBJECT_H
 
-// STD includes
+// STL includes
 #include <vector>
 #include <memory>
+#include <unordered_map>
 
 // SimMedTK includes
 #include "Core/Vector.h"
@@ -122,14 +123,102 @@ public:
     /// \brief print information related the scene object
     virtual void printInfo() const = 0;
 
+    ///
+    /// \brief Set to compute contact forces
+    ///
+    bool computeContactForce()
+    {
+        return this->hasContactForces;
+    }
+
+    ///
+    /// \brief Set to not compute contact forces
+    ///
+    void setContactForcesOff()
+    {
+        this->hasContactForces = false;
+    }
+
+    ///
+    /// \brief Set to not compute contact forces
+    ///
+    void setContactForcesOn()
+    {
+        this->hasContactForces = true;
+    }
+
+    // Get contact forces vector
+    std::unordered_map<int,core::Vec3d> &getContactForces()
+    {
+        return this->contactForces;
+    }
+
+    const std::unordered_map<int,core::Vec3d> &getContactForces() const
+    {
+        return this->contactForces;
+    }
+
+    // Get contact forces vector
+    std::unordered_map<int,core::Vec3d> &getContactPoints()
+    {
+        return this->contactPoints;
+    }
+
+    const std::unordered_map<int,core::Vec3d> &getContactPoints() const
+    {
+        return this->contactPoints;
+    }
+
+    /// \brief  returns velocity of at a given location
+    /// (not given node) in contact force vector
+    virtual core::Vec3d getVelocity(const int) const
+    {
+        return core::Vec3d::Zero();
+    }
+
+    /// \brief Set all contact forces to zero (if any)
+    void setContactForcesToZero()
+    {
+        this->contactForces.clear();
+    }
+
+    void setContactForce(const int dofID,const core::Vec3d &force)
+    {
+        this->contactForces[dofID] = force;
+    }
+
+    void setContactForce(const int dofID,
+                         const core::Vec3d &point,
+                         const core::Vec3d &force)
+    {
+        this->contactPoints[dofID] = point;
+        this->contactForces[dofID] = force;
+    }
+
+
+    void setModel(std::shared_ptr<Model> m)
+    {
+        this->model = m;
+    }
+
+    std::shared_ptr<Model> getModel()
+    {
+        return this->model;
+    }
+
 protected:
     bool isActive;
 
 private:
-    std::shared_ptr<ObjectSimulator> objectSim; // object simulator that will simulate the object
+    std::shared_ptr<ObjectSimulator> objectSim; //!< object simulator that will simulate the object
     std::shared_ptr<CustomRenderer> customRender;
-    std::vector<core::Vec3d> localVertices; // local copy of vertices
+    std::vector<core::Vec3d> localVertices; //!< local copy of vertices
     ObjectInitFlags flags;
+    bool hasContactForces;
+    std::unordered_map<int,core::Vec3d> contactForces;
+    std::unordered_map<int,core::Vec3d> contactPoints;
+
+    std::shared_ptr<Model> model; //!< model attached to this scene object
 };
 
 #endif
diff --git a/Core/StaticSceneObject.cpp b/Core/StaticSceneObject.cpp
index 7d6bbbbc8996ea3b252306d170bd9c8cabf9f8a0..19724c22ee49112eb857fac2bc65ff4350d7b115 100644
--- a/Core/StaticSceneObject.cpp
+++ b/Core/StaticSceneObject.cpp
@@ -66,16 +66,6 @@ std::shared_ptr<SceneObject> StaticSceneObject::clone()
     return safeDownCast<SceneObject>();
 }
 
-void StaticSceneObject::setModel(std::shared_ptr<Model> model)
-{
-    this->staticModel = model;
-}
-
-std::shared_ptr<Model> StaticSceneObject::getModel()
-{
-    return staticModel;
-}
-
 void StaticSceneObject::printInfo() const
 {
     std::cout << "\t-------------------------------------\n";
diff --git a/Core/StaticSceneObject.h b/Core/StaticSceneObject.h
index 283872761f51c7579ccbb4f9f89ba94537d5577c..eab0587609ae6b27362baa324f762d8f21089945 100644
--- a/Core/StaticSceneObject.h
+++ b/Core/StaticSceneObject.h
@@ -24,9 +24,6 @@
 #ifndef SMSTATICSCENEOBJECT_H
 #define SMSTATICSCENEOBJECT_H
 
-// STL includes
-#include <unordered_map>
-
 // SimMedTK includes
 #include "Core/Config.h"
 #include "Core/Model.h"
@@ -73,37 +70,6 @@ public:
     void printInfo() const override;
 
     virtual void handleEvent(std::shared_ptr<core::Event>) override {}
-
-    void setModel(std::shared_ptr<Model> model);
-
-    std::shared_ptr<Model> getModel();
-
-    // Get contact forces vector
-    std::unordered_map<int,core::Vec3d> &getContactForces()
-    {
-        return this->contactForces;
-    }
-
-    const std::unordered_map<int,core::Vec3d> &getContactForces() const
-    {
-        return this->contactForces;
-    }
-
-    // Get contact forces vector
-    std::unordered_map<int,core::Vec3d> &getContactPoints()
-    {
-        return this->contactPoints;
-    }
-
-    const std::unordered_map<int,core::Vec3d> &getContactPoints() const
-    {
-        return this->contactPoints;
-    }
-public:
-    /// \brief static scene object contains a mesh
-    std::shared_ptr<Model> staticModel;
-    std::unordered_map<int,core::Vec3d> contactForces;
-    std::unordered_map<int,core::Vec3d> contactPoints;
 };
 
 #endif
diff --git a/Examples/vegaFem/CMakeLists.txt b/Examples/vegaFem/CMakeLists.txt
index ac36b6a40d6b05accbf6b22a0b341bf92b18c168..0c54649d67cb0b5afa81320d1a57c4d8a7cad179 100644
--- a/Examples/vegaFem/CMakeLists.txt
+++ b/Examples/vegaFem/CMakeLists.txt
@@ -13,6 +13,8 @@ target_link_libraries(${APP}
   ContactHandling
   Event
   IO
+  Devices
+  VirtualTools
 )
 
 set_target_properties(${APP}
diff --git a/Examples/vegaFem/main.cpp b/Examples/vegaFem/main.cpp
index e1981ccb84762df30ad61f71448d58a8a009ea79..a29edf6f7727930b59dd15aa077e312d76b37563 100644
--- a/Examples/vegaFem/main.cpp
+++ b/Examples/vegaFem/main.cpp
@@ -67,19 +67,19 @@ int main(int ac, char** av)
     // 3. Create default scene (scene 0)
     //-------------------------------------------------------
     auto sdk = SDK::createStandardSDK();
-    auto client = std::make_shared<VRPNForceDevice>();
-    //get some user input and setup device url
-    std::string input = "navigator@localhost";
-    std::cout << "Enter the VRPN device URL(" << client->getDeviceURL() << "): ";
-    std::getline(std::cin, input);
-    if(!input.empty())
-    {
-        client->setDeviceURL(input);
-    }
-    auto controller = std::make_shared<ToolCoupler>(client);
-    controller->setScalingFactor(5.0);
-    sdk->registerModule(client);
-    sdk->registerModule(controller);
+//     auto client = std::make_shared<VRPNForceDevice>();
+//     //get some user input and setup device url
+//     std::string input = "navigator@localhost";
+//     std::cout << "Enter the VRPN device URL(" << client->getDeviceURL() << "): ";
+//     std::getline(std::cin, input);
+//     if(!input.empty())
+//     {
+//         client->setDeviceURL(input);
+//     }
+//     auto controller = std::make_shared<ToolCoupler>(client);
+//     controller->setScalingFactor(5.0);
+//     sdk->registerModule(client);
+//     sdk->registerModule(controller);
 
     //-------------------------------------------------------
     // Create scene actor 1:  fem scene object + fem simulator
@@ -87,11 +87,12 @@ int main(int ac, char** av)
 
     // create a FEM simulator
     auto femSimulator = std::make_shared<VegaFemSimulator>(sdk->getErrorLog());
-    femSimulator->setHapticTool(controller);
+//     femSimulator->setHapticTool(controller);
 
     // create a Vega based FEM object and attach it to the fem simulator
     auto femObject = std::make_shared<VegaFemSceneObject>(
         sdk->getErrorLog(),configFile);
+    femObject->setContactForcesOn();
 
     auto meshRenderDetail = std::make_shared<RenderDetail>(SIMMEDTK_RENDER_WIREFRAME |
     //| SIMMEDTK_RENDER_VERTICES
@@ -131,21 +132,21 @@ int main(int ac, char** av)
     // Create scene actor 2:  loli tool
     // create a static object to hold the lolitool scene object of given normal and position
     //-------------------------------------------------------
-    auto loliSceneObject = std::make_shared<StaticSceneObject>();
-
-    auto loliCollisionModel = std::make_shared<MeshCollisionModel>();
-    loliCollisionModel->loadTriangleMesh("./loli.vtk");
-    loliSceneObject->setModel(loliCollisionModel);
-
-    auto loliMesh = loliCollisionModel->getMesh();
-    Core::BaseMesh::TransformType transform = Eigen::Translation3d(core::Vec3d(0,5,0))*Eigen::Scaling(.5);
-
-    loliMesh->transform(transform);
-    loliMesh->updateInitialVertices();
-
-    auto loliSimulator = std::make_shared<DefaultSimulator>(sdk->getErrorLog());
-    sdk->addSceneActor(loliSceneObject, loliSimulator);
-    controller->setMesh(loliCollisionModel->getMesh());
+//     auto loliSceneObject = std::make_shared<StaticSceneObject>();
+//
+//     auto loliCollisionModel = std::make_shared<MeshCollisionModel>();
+//     loliCollisionModel->loadTriangleMesh("./loli.vtk");
+//     loliSceneObject->setModel(loliCollisionModel);
+//
+//     auto loliMesh = loliCollisionModel->getMesh();
+//     Core::BaseMesh::TransformType transform = Eigen::Translation3d(core::Vec3d(0,5,0))*Eigen::Scaling(.5);
+//
+//     loliMesh->transform(transform);
+//     loliMesh->updateInitialVertices();
+//
+//     auto loliSimulator = std::make_shared<DefaultSimulator>(sdk->getErrorLog());
+//     sdk->addSceneActor(loliSceneObject, loliSimulator);
+//     controller->setMesh(loliCollisionModel->getMesh());
 
     //-------------------------------------------------------
     // Register both object simulators
diff --git a/Simulators/SceneObjectDeformable.cpp b/Simulators/SceneObjectDeformable.cpp
index b706630fd974d14ca237892abd391d884ccbc202..be957d55b9c575c28cf66e794d40b65076757d56 100644
--- a/Simulators/SceneObjectDeformable.cpp
+++ b/Simulators/SceneObjectDeformable.cpp
@@ -45,7 +45,7 @@ SceneObjectDeformable::~SceneObjectDeformable()
 
 void SceneObjectDeformable::applyContactForces()
 {
-    for(const auto &cf : contactForces)
+    for(const auto &cf : this->getContactForces())
     {
         int i = cf.first;
         f_ext[i] += cf.second[0];
@@ -54,27 +54,8 @@ void SceneObjectDeformable::applyContactForces()
     }
 }
 
-void SceneObjectDeformable::setContactForcesToZero()
-{
-    contactForces.clear();
-    contactPoints.clear();
-}
-
-void SceneObjectDeformable::setContactForceOfNodeWithDofID(const int dofID,
-                                                             const core::Vec3d &force)
-{
-    contactForces[dofID] = force;
-}
-
-void SceneObjectDeformable::setContactForceOfNodeWithDofID(const int dofID,
-                                                           const core::Vec3d &point,
-                                                           const core::Vec3d &force)
-{
-    contactPoints[dofID] = point;
-    contactForces[dofID] = force;
-}
 
-core::Vec3d SceneObjectDeformable::getVelocityOfNodeWithDofID(const int dofID) const
+core::Vec3d SceneObjectDeformable::getVelocity(const int dofID) const
 {
     core::Vec3d vel(uvel[dofID], uvel[dofID + 1], uvel[dofID + 2]);
 
diff --git a/Simulators/SceneObjectDeformable.h b/Simulators/SceneObjectDeformable.h
index 16f244f8f2d659e9fddab739358c75e64929debe..2a21bdbed4e7914b01abafc503b9a5b7905b99bf 100644
--- a/Simulators/SceneObjectDeformable.h
+++ b/Simulators/SceneObjectDeformable.h
@@ -72,25 +72,13 @@ public:
     /// \brief Append the contact forces (if any) to external forces
     void applyContactForces();
 
-    /// \brief Set all contact forces to zero (if any)
-    void setContactForcesToZero();
-
-    /// \brief  Sets the contact force at a given location
-    /// (not given node) in contact force vector
-    void setContactForceOfNodeWithDofID(const int dofID, const core::Vec3d &force);
-
-    /// \brief  Sets the contact force at a given location
-    /// (not given node) in contact force vector
-    void setContactForceOfNodeWithDofID(const int dofID, const core::Vec3d &contactPoint,
-                                        const core::Vec3d &force);
-
     /// \brief  returns displacement of at a given location
     /// (not given node) in contact force vector
     core::Vec3d getDisplacementOfNodeWithDofID(const int dofID) const;
 
     /// \brief  returns velocity of at a given location
     /// (not given node) in contact force vector
-    core::Vec3d getVelocityOfNodeWithDofID(const int dofID) const;
+    core::Vec3d getVelocity(const int dofID) const;
 
     /// \brief  returns acceleration of at a given location
     /// (not given node) in contact force vector
@@ -150,28 +138,6 @@ public:
         return this->f_ext;
     }
 
-    // Get contact forces vector
-    std::unordered_map<int,core::Vec3d> &getContactForces()
-    {
-        return this->contactForces;
-    }
-
-    const std::unordered_map<int,core::Vec3d> &getContactForces() const
-    {
-        return this->contactForces;
-    }
-
-    // Get contact forces vector
-    std::unordered_map<int,core::Vec3d> &getContactPoints()
-    {
-        return this->contactPoints;
-    }
-
-    const std::unordered_map<int,core::Vec3d> &getContactPoints() const
-    {
-        return this->contactPoints;
-    }
-
 protected:
     friend class SceneObjectDeformableRenderDelegate;
 
@@ -205,8 +171,6 @@ protected:
     std::shared_ptr<SurfaceMesh> primarySurfaceMesh;
     std::shared_ptr<SurfaceMesh> secondarySurfaceMesh;
 
-    std::unordered_map<int,core::Vec3d> contactForces;
-    std::unordered_map<int,core::Vec3d> contactPoints;
 };
 
 #endif // SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
diff --git a/Simulators/VegaFemSceneObject.cpp b/Simulators/VegaFemSceneObject.cpp
index 5811ebca915123f14d149ed5d517882c2a56677b..452e832a6988f461fbcb9c2fa8da3275e979b7ea 100644
--- a/Simulators/VegaFemSceneObject.cpp
+++ b/Simulators/VegaFemSceneObject.cpp
@@ -1002,7 +1002,7 @@ void VegaFemSceneObject::printInfo() const
     std::cout << "\t-------------------------------------\n";
 }
 
-core::Vec3d VegaFemSceneObject::getVelocityOfNodeWithDofID(const int dofID) const
+core::Vec3d VegaFemSceneObject::getVelocity(const int dofID) const
 {
     core::Vec3d vel(uvel[dofID], uvel[dofID + 1], uvel[dofID + 2]);
 
diff --git a/Simulators/VegaFemSceneObject.h b/Simulators/VegaFemSceneObject.h
index 8bb0a65132f5fe8c971a2ff1026f1ca05305db70..9d4c6f6b940c8b79c75657cf8e058209f49dd434 100644
--- a/Simulators/VegaFemSceneObject.h
+++ b/Simulators/VegaFemSceneObject.h
@@ -152,7 +152,7 @@ public:
 
     /// \brief returns velocity given the
     /// localtion in the global velocity vector
-    core::Vec3d getVelocityOfNodeWithDofID(const int dofID) const;
+    core::Vec3d getVelocity(const int dofID) const;
 
     /// \brief returns displacement given the
     /// localtion in the global displacement vector
diff --git a/VirtualTools/CMakeLists.txt b/VirtualTools/CMakeLists.txt
index d8dc747daf60fd75a5854dbb46a50d840b605a6b..ba2c38d45b896f22aa6a41523c0386b8933381a9 100644
--- a/VirtualTools/CMakeLists.txt
+++ b/VirtualTools/CMakeLists.txt
@@ -12,7 +12,6 @@ target_link_libraries(VirtualTools
   PRIVATE
     Core
     Mesh
-    Simulators
     Event
     Devices
 )