diff --git a/src/Tearing/BaseTearingEngine.h b/src/Tearing/BaseTearingEngine.h
index c853a8d..9d9b8e5 100644
--- a/src/Tearing/BaseTearingEngine.h
+++ b/src/Tearing/BaseTearingEngine.h
@@ -125,16 +125,16 @@ class BaseTearingEngine : public core::DataEngine
 
 	virtual void computeFracturePath() = 0;
 	
-	void computeFractureDirection(const Coord principleStressDirection, Coord& fracture_direction);
+	Coord computeFractureDirection(const Coord& principleStressDirection);
 
 	/// <summary>
 	/// compute extremities of fracture Pb and Pc from a start point Pa
 	/// </summary>
 	/// @param Pa - point with maxStress where fracture start
-	/// @param direction - direction of maximum principal stress
+	/// @param fractureDirection - direction of fracture
 	/// @return Pb - one of the extremities of fracture
 	/// @return Pc - one of the extremities of fracture
-	virtual void computeEndPoints(Coord Pa, Coord direction, Coord& Pb, Coord& Pc);
+	virtual void computeEndPoints(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc);
 
 	/// <summary>
 	/// compute ignored triangle at start of the tearing algo
diff --git a/src/Tearing/BaseTearingEngine.inl b/src/Tearing/BaseTearingEngine.inl
index 5925916..4778232 100644
--- a/src/Tearing/BaseTearingEngine.inl
+++ b/src/Tearing/BaseTearingEngine.inl
@@ -329,19 +329,20 @@ void BaseTearingEngine<DataTypes>::updateTriangleInformation()
 
 
 template<class DataTypes>
-inline void BaseTearingEngine<DataTypes>::computeFractureDirection(const Coord principleStressDirection,Coord & fracture_direction)
+typename DataTypes::Coord BaseTearingEngine<DataTypes>::computeFractureDirection(const Coord& principleStressDirection)
 {
+    Coord fracture_direction = { 0.0, 0.0, 0.0 };
+
     if (m_maxStressTriangleIndex == InvalidID) {
-        fracture_direction = { 0.0, 0.0, 0.0 };
-        return;
+        return fracture_direction;
     }
 
     const Triangle& VertexIndicies = m_topology->getTriangle(m_maxStressTriangleIndex);
-    constexpr size_t numVertices = 3;
 
-    Index B_id = -1, C_id = -1;
+    Index B_id = sofa::InvalidID;
+    Index C_id = sofa::InvalidID;
 
-    for (unsigned int vertex_id = 0; vertex_id < numVertices; vertex_id++)
+    for (sofa::Index vertex_id = 0; vertex_id < 3; vertex_id++)
     {
         if (VertexIndicies[vertex_id] == m_maxStressVertexIndex)
         {
@@ -356,27 +357,22 @@ inline void BaseTearingEngine<DataTypes>::computeFractureDirection(const Coord p
     Coord B = x[B_id];
     Coord C = x[C_id];
 
-    Coord AB = B - A;
-    Coord AC = C - A;
-
-    Coord triangleNormal = sofa::type::cross(AB,AC);
+    Coord triangleNormal = sofa::type::cross(B - A, C - A);
     fracture_direction = sofa::type::cross(triangleNormal, principleStressDirection);
+
+    return fracture_direction;
 }
 
 
 template <class DataTypes>
-void BaseTearingEngine<DataTypes>::computeEndPoints(
-    Coord Pa,
-    Coord direction,
-    Coord& Pb, Coord& Pc)
+void BaseTearingEngine<DataTypes>::computeEndPoints(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc)
 {
-    Coord fractureDirection;
-    computeFractureDirection(direction, fractureDirection);
     Real norm_fractureDirection = fractureDirection.norm();
     Pb = Pa + d_fractureMaxLength.getValue() / norm_fractureDirection * fractureDirection;
     Pc = Pa - d_fractureMaxLength.getValue() / norm_fractureDirection * fractureDirection;
 }
 
+
 template <class DataTypes>
 void BaseTearingEngine<DataTypes>::computeTriangleToSkip()
 {
diff --git a/src/Tearing/TearingAlgorithms.h b/src/Tearing/TearingAlgorithms.h
index 5bd4157..da4037b 100644
--- a/src/Tearing/TearingAlgorithms.h
+++ b/src/Tearing/TearingAlgorithms.h
@@ -74,6 +74,16 @@ class TearingAlgorithms
 	
 	void computeFracturePath(FracturePath& my_fracturePath);
 
+
+	/// <summary>
+	/// computes the extremities of the (normalized) fracture PbPa on the edge of the triangle
+	/// </summary>
+	/// @param Pa - the point where the fracture starts
+	/// @param normalizedFractureDirection - normalized fracture direction
+	/// @return Pb - one of the extremities of fracture
+	/// @return t - a parameter needed to calculate Pb
+	TriangleID computeIntersectionNeighborTriangle(const Index ptAId, const Coord& ptA, const Coord& normalizedFractureDirection, Coord& Pb);
+
 	void computeFracturePath(const Coord& pA, Index triId, const Coord pB, const Coord pC);
 
 	/// <summary>
diff --git a/src/Tearing/TearingAlgorithms.inl b/src/Tearing/TearingAlgorithms.inl
index c4fd82e..c010022 100644
--- a/src/Tearing/TearingAlgorithms.inl
+++ b/src/Tearing/TearingAlgorithms.inl
@@ -51,6 +51,67 @@ TearingAlgorithms<DataTypes>::~TearingAlgorithms()
 }
 
 
+template<class DataTypes>
+inline TriangleID TearingAlgorithms<DataTypes>::computeIntersectionNeighborTriangle(const Index ptAId, const Coord& ptA, const Coord& normalizedFractureDirection, Coord& ptB)
+{
+    // Get triangle in the fracture direction
+    TriangleID theTriId = m_triangleGeo->getTriangleInDirection(ptAId, normalizedFractureDirection);
+    
+    // If not, could be on the border. Return invalidID
+    if (theTriId > this->m_topology->getNbTriangles() - 1)
+        return sofa::InvalidID;
+
+    // If it exists, get triangle and search for ptAId local index in the triangle
+    const Triangle& theTri = this->m_topology->getTriangle(theTriId);
+    PointID localAId = sofa::InvalidID;
+    for (PointID vertex_id = 0; vertex_id < 3; vertex_id++)
+    {
+        if (theTri[vertex_id] == ptAId)
+        {
+            localAId = vertex_id;
+            break;
+        }
+    }
+
+    // Get the opposite edge
+    EdgeID oppositeEdgeId = this->m_topology->getEdgesInTriangle(theTriId)[localAId];
+    Edge oppositeEdge = this->m_topology->getEdge(oppositeEdgeId);
+
+
+    //Building point B such that to be sure that AB intersects CD, based on "Losange"
+    const Coord pE0 = this->m_triangleGeo->getPointPosition(oppositeEdge[0]);
+    const Coord pE1 = this->m_triangleGeo->getPointPosition(oppositeEdge[1]);
+
+    const Real AC_length = (pE0 - ptA).norm();
+    const Real AD_length = (pE1 - ptA).norm();
+    const Real Length = AC_length + AD_length;
+
+    ptB = ptA + Length * normalizedFractureDirection;
+
+    // Compute intersection on the opposite edge
+    sofa::type::vector<EdgeID> intersectedEdges;
+    sofa::type::vector<Real> baryCoefs;
+    m_triangleGeo->computeSegmentTriangleIntersectionInPlane(ptA, ptB, theTriId, intersectedEdges, baryCoefs);
+
+    bool found = false;
+    for (unsigned int i=0; i< intersectedEdges.size(); ++i)
+    {
+        if (intersectedEdges[i] == oppositeEdgeId)
+        {
+            found = true;
+            ptB = pE0 * baryCoefs[i] + pE1 * (1 - baryCoefs[i]);
+            break;
+        }
+    }
+
+    if (!found)
+        return sofa::InvalidID;
+    else
+        return theTriId;
+}
+
+
+
 template <class DataTypes>
 void TearingAlgorithms<DataTypes>::computeFracturePath(const Coord& pA, Index triId, const Coord pB, const Coord pC)
 {
diff --git a/src/Tearing/TearingEngine.h b/src/Tearing/TearingEngine.h
index d56342d..6436e8c 100644
--- a/src/Tearing/TearingEngine.h
+++ b/src/Tearing/TearingEngine.h
@@ -81,20 +81,11 @@ class TearingEngine : public BaseTearingEngine<DataTypes>
 	/// computes the extremities of fracture Pb and Pc on the edge of neighboring triangles
 	/// </summary>
 	/// @param Pa - the point where the fracture starts
-	/// @param direction - principle stress direction
+	/// @param fractureDirection - fracture direction
 	/// @return Pb - one of the extremities of fracture
 	/// @return Pc - one of the extremities of fracture
-	bool computeEndPointsNeighboringTriangles(Coord Pa, Coord direction, Coord& Pb, Coord& Pc);
+	bool computeEndPointsNeighboringTriangles(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc);
 
-	/// <summary>
-	/// computes the extremities of the (normalized) fracture PbPa on the edge of the triangle
-	/// </summary>
-	/// @param Pa - the point where the fracture starts
-	/// @param normalizedFractureDirection - normalized fracture direction
-	/// @return Pb - one of the extremities of fracture
-	/// @return t - a parameter needed to calculate Pb
-	bool computeIntersectionNeighborTriangle(Coord normalizedFractureDirection, Coord Pa, Coord& Pb, Real& t);
-	
 	void algoFracturePath() override;
 
 	void computeFracturePath() override;
diff --git a/src/Tearing/TearingEngine.inl b/src/Tearing/TearingEngine.inl
index 5206ac6..6e9585c 100644
--- a/src/Tearing/TearingEngine.inl
+++ b/src/Tearing/TearingEngine.inl
@@ -41,85 +41,26 @@ using sofa::type::Vec3;
 // --------------------------------------------------------------------------------------
 // --- Computation methods
 // --------------------------------------------------------------------------------------
-template<class DataTypes>
-inline bool TearingEngine<DataTypes>::computeIntersectionNeighborTriangle(Coord normalizedFractureDirection, Coord Pa, Coord& Pb, Real& t)
-{
-    SOFA_UNUSED(Pa);
-
-    if (m_maxStressVertexIndex == InvalidID)
-        return false;    
-
-    // Get Geometry Algorithm
-    TriangleSetGeometryAlgorithms<DataTypes>* _triangleGeo = nullptr;
-    this->m_topology->getContext()->get(_triangleGeo);
-    if (!_triangleGeo)
-    {
-        msg_error() << "Missing component: Unable to get TriangleSetGeometryAlgorithms from the current context.";
-        sofa::core::objectmodel::BaseObject::d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
-        return false;
-    }
-
-
-    Index triangle_id = _triangleGeo->getTriangleInDirection(m_maxStressVertexIndex, normalizedFractureDirection);
-    if (triangle_id > this->m_topology->getNbTriangles() - 1)
-        return false;
-
-
-    const Triangle& VertexIndicies = this->m_topology->getTriangle(triangle_id);
-
-    constexpr size_t numVertices = 3;
-    Index B_id = -1, C_id = -1;
-
-    for (unsigned int vertex_id = 0; vertex_id < numVertices; vertex_id++)
-    {
-        if (VertexIndicies[vertex_id] == m_maxStressVertexIndex)
-        {
-            B_id = VertexIndicies[(vertex_id + 1) % 3];
-            C_id = VertexIndicies[(vertex_id + 2) % 3];
-            break;
-        }
-
-    }
-
-    helper::ReadAccessor< Data<VecCoord> > x(d_input_positions);
-    Coord A = x[m_maxStressVertexIndex];
-    Coord B = x[B_id];
-    Coord C = x[C_id];
-
-    if (rayTriangleIntersection(A, B, C, normalizedFractureDirection, t, Pb))
-        return true;
-    else
-        return false;
-
-}
 
 template<class DataTypes>
-inline bool TearingEngine<DataTypes>::computeEndPointsNeighboringTriangles(Coord Pa, Coord direction, Coord& Pb, Coord& Pc)
+inline bool TearingEngine<DataTypes>::computeEndPointsNeighboringTriangles(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc)
 {
     bool t_b_ok = false; 
     bool t_c_ok = false;
-    //compute fracture direction perpendicular to the principal stress direction
-    Coord fractureDirection;
-    this->computeFractureDirection(direction, fractureDirection);
    
-    Real norm_fractureDirection = fractureDirection.norm();
-    Coord dir_b = 1.0 / norm_fractureDirection * fractureDirection;
-
-    Real t_b;
-    if (computeIntersectionNeighborTriangle(dir_b,Pa, Pb, t_b))
-        t_b_ok = true;
+    Coord dir_b = fractureDirection / fractureDirection.norm();
+    TriangleID t_b = this->m_tearingAlgo->computeIntersectionNeighborTriangle(m_maxStressVertexIndex, Pa, dir_b, Pb);
 
-    
     Coord dir_c = -dir_b;
-    Real t_c;
-    if (computeIntersectionNeighborTriangle(dir_c,Pa,Pc, t_c))
-        t_c_ok = true;
+    TriangleID t_c = this->m_tearingAlgo->computeIntersectionNeighborTriangle(m_maxStressVertexIndex, Pa, dir_c, Pc);
+
 
-    if (!(t_b_ok && t_c_ok))
+    if (t_b == sofa::InvalidID || t_c == sofa::InvalidID)
     {
         msg_warning() << "Not able to build the fracture path through neighboring triangles.";
         return false;
     }
+
     return true;
 }
 
@@ -163,33 +104,36 @@ void TearingEngine<DataTypes>::computeFracturePath()
     // frist clear ereything
     this->clearFracturePath();
 
-    if (m_maxStressTriangleIndex != InvalidID) // we have triangle to start and also a vertex id
+    if (m_maxStressTriangleIndex == InvalidID) // no candidate to start fracture
+        return;
+
+    //Recording the endpoints of the fracture segment
+    helper::ReadAccessor< Data<VecCoord> > x(d_input_positions);
+
+    const Coord fractureDirection = this->computeFractureDirection(m_triangleInfoTearing[m_maxStressTriangleIndex].principalStressDirection);
+
+    const Coord Pa = x[m_maxStressVertexIndex];
+    Coord Pb, Pc;
+
+    if (this->d_fractureMaxLength.getValue() == 0.0) 
+    {
+        bool checkEndsPoints = this->computeEndPointsNeighboringTriangles(Pa, fractureDirection, Pb, Pc);
+        if (!checkEndsPoints)
+            return;
+    }
+    else
     {
-        //Recording the endpoints of the fracture segment
-        helper::ReadAccessor< Data<VecCoord> > x(d_input_positions);
-
-        Coord principalStressDirection = m_triangleInfoTearing[m_maxStressTriangleIndex].principalStressDirection;
-        Coord Pa = x[m_maxStressVertexIndex];      
-        Coord Pb, Pc;
-
-        if (this->d_fractureMaxLength.getValue() == 0.0) {
-            computeEndPointsNeighboringTriangles(Pa, principalStressDirection, Pb, Pc);
-        }
-        else 
-        {
-            this->computeEndPoints(Pa, principalStressDirection, Pb, Pc);
-        }
-
-        this->m_fracturePath.ptA = type::Vec3(Pa[0], Pa[1], Pa[2]);
-        this->m_fracturePath.ptB = type::Vec3(Pb[0], Pb[1], Pb[2]);
-        this->m_fracturePath.ptC = type::Vec3(Pc[0], Pc[1], Pc[2]);
-        this->m_fracturePath.triIdA = m_maxStressTriangleIndex;
-
-        this->m_tearingAlgo->computeFracturePath(this->m_fracturePath);
-        this->m_stepCounter++;
-
-        //this->m_tearingAlgo->computeFracturePath(Pa, m_maxStressTriangleIndex, Pb, Pc);
+        this->computeEndPoints(Pa, fractureDirection, Pb, Pc); // compute orthogonal fracture using d_fractureMaxLength
     }
+    
+
+    this->m_fracturePath.ptA = type::Vec3(DataTypes::getCPos(Pa));
+    this->m_fracturePath.ptB = type::Vec3(DataTypes::getCPos(Pb));
+    this->m_fracturePath.ptC = type::Vec3(DataTypes::getCPos(Pc));
+    this->m_fracturePath.triIdA = m_maxStressTriangleIndex;
+
+    this->m_tearingAlgo->computeFracturePath(this->m_fracturePath);
+    this->m_stepCounter++;
 }
 
 
diff --git a/src/Tearing/TearingScenarioEngine.h b/src/Tearing/TearingScenarioEngine.h
index 59c2d4f..3913d77 100644
--- a/src/Tearing/TearingScenarioEngine.h
+++ b/src/Tearing/TearingScenarioEngine.h
@@ -77,7 +77,7 @@ class TearingScenarioEngine : public BaseTearingEngine<DataTypes>
 	void algoFracturePath() override;
 	void computeFracturePath() override {};
 
-	void computeEndPoints(Coord Pa, Coord direction, Coord& Pb, Coord& Pc) override;
+	void computeEndPoints(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc) override;
 
 	/// Value to store scenario fracture path
 	Coord m_Pa, m_Pb, m_Pc;
diff --git a/src/Tearing/TearingScenarioEngine.inl b/src/Tearing/TearingScenarioEngine.inl
index 18fce34..0dd17a8 100644
--- a/src/Tearing/TearingScenarioEngine.inl
+++ b/src/Tearing/TearingScenarioEngine.inl
@@ -106,17 +106,14 @@ inline void TearingScenarioEngine<DataTypes>::computePerpendicular(Coord dir, Co
 }
 
 template <class DataTypes>
-void TearingScenarioEngine<DataTypes>::computeEndPoints(
-    Coord Pa,
-    Coord dir,
-    Coord& Pb, Coord& Pc)
+void TearingScenarioEngine<DataTypes>::computeEndPoints(const Coord& Pa, const Coord& fractureDirection, Coord& Pb, Coord& Pc)
 {
     const Real& alpha = d_startLength.getValue();
         
-    Real norm_dir = dir.norm();
+    Real norm_dir = fractureDirection.norm();
        
-    Pb = Pa + alpha/norm_dir * dir;
-    Pc = Pa - alpha /norm_dir * dir;
+    Pb = Pa + alpha/norm_dir * fractureDirection;
+    Pc = Pa - alpha /norm_dir * fractureDirection;
 }
 
 template <class DataTypes>