[src] Fix the logic when to start tearing and add fracture structure (#65)

* [src] Fix the logic when to start tearing

* Backup work on new fracture path

* clean PR

Author: epernod

scenes/TearingEngine/SquareTissue_11_MiddlePressure .scn

@@ -39,9 +39,9 @@
         <CGLinearSolver iterations="25" name="linear solver" tolerance="1.0e-9" threshold="1.0e-9" />
         <MeshGmshLoader name="meshLoader" filename="mesh/square3.msh" scale="10" createSubelements="true" />
 
-		<TriangleSetTopologyContainer  name="Container" src="@meshLoader" fileTopology="" tags=" " />
+		<TriangleSetTopologyContainer  name="Container" src="@meshLoader"  />
 		<TriangleSetTopologyModifier   name="Modifier" />
-		<TriangleSetGeometryAlgorithms name="GeomAlgo"   template="Vec3d" />
+		<TriangleSetGeometryAlgorithms name="GeomAlgo"   template="Vec3d" showEdgeIndices="1" showTriangleIndices="1" showIndicesScale="0.001"/>
 
         <MechanicalObject name="Mo"/>
         <DiagonalMass massDensity="0.15" />
@@ -51,15 +51,15 @@
 
 		<BoxROI template="Vec3d" name="PullMiddle" box="4.9 4.95 0.1  5.1 5.1 0.1" drawBoxes="1"/>
 		<!--<BoxROI template="Vec3d" name="PulledSide" box="-0.1 0.1 0.1  10.1 -0.1 -0.1" drawBoxes="1"/>-->
-		<ConstantForceField name="CFF" indices="772" forces="0 -300 0" showArrowSize=".01" />        
+		<ConstantForceField name="CFF" indices="772" forces="0 -500 0" showArrowSize=".01" />        
 
 		<TriangularFEMForceField name="FEM" youngModulus="300" poissonRatio="0.3" method="large" computePrincipalStress="true"/>
         <TriangularBendingSprings name="FEM-Bend" stiffness="300" damping="1.0" />
         <PlaneForceField normal="0 1 0" d="-30" stiffness="100000" showPlane="1" showPlaneSize="20"/>
 
         <TriangleCollisionModel />
-		<TearingEngine name="TearingEngine" input_position="@Mo.position" stressThreshold="0.1" 
-            step="20" nbFractureMax="40" showFracturePath="1"/>
+		<TearingEngine name="TearingEngine" input_position="@Mo.position" stressThreshold="15.0" 
+            step="20" nbFractureMax="40" showFracturePath="1" showTearableCandidates="1" fractureMaxLength="1.0"/>
 
         <Node >
             <OglModel name="Visual" color="gray" />

scenes/TearingEngine/SquareTissue_11_MiddlePressure .scn.qglviewer.view

@@ -1,2 +1,2 @@
-2.39377 7.2768 19.272
+5.18732 3.91825 1.3
 -0.0438528 -0.046858 0.00673547 0.997916

src/Tearing/BaseTearingEngine.inl

@@ -162,7 +162,7 @@ void BaseTearingEngine<DataTypes>::doUpdate()
     if (sofa::core::objectmodel::BaseObject::d_componentState.getValue() != sofa::core::objectmodel::ComponentState::Valid)
         return;
 
-    m_stepCounter++;   
+    
 
     if (d_ignoreTriangles.getValue())
     {
@@ -215,7 +215,6 @@ void BaseTearingEngine<DataTypes>::triangleOverThresholdPrincipalStress()
     candidate.clear();
     maxStress = 0;
     helper::WriteAccessor< Data<vector<Index>> >triangleToSkip(d_trianglesToIgnore);
-    
     m_maxStressTriangleIndex = InvalidID;
     for (unsigned int i = 0; i < triangleList.size(); i++)
     {
@@ -639,6 +638,9 @@ void BaseTearingEngine<DataTypes>::handleEvent(sofa::core::objectmodel::Event* e
         return; // We only launch computation at end of a simulation step
     }
 
+    if (m_tearingAlgo->getFractureNumber() > d_nbFractureMax.getValue())
+        return;
+
     computeFracturePath();
 
     // Hack: we access one output value to force the engine to call doUpdate()
@@ -647,14 +649,9 @@ void BaseTearingEngine<DataTypes>::handleEvent(sofa::core::objectmodel::Event* e
 
     // Perform fracture every d_stepModulo
     int step = d_stepModulo.getValue();
-    if (step == 0) // interactive version
-    {
-        if (m_stepCounter > 200 && (m_tearingAlgo->getFractureNumber() < d_nbFractureMax.getValue())){
-            algoFracturePath();
-        }
-    }
-    else if (((m_stepCounter % step) == 0) && (m_tearingAlgo->getFractureNumber() < d_nbFractureMax.getValue()))
+    if (m_stepCounter > step)
     {
+        m_stepCounter = 0;
         algoFracturePath();
     }
 }

src/Tearing/TearingAlgorithms.h

@@ -29,6 +29,7 @@
 #include <sofa/component/topology/container/dynamic/TriangleSetTopologyModifier.h>
 #include <sofa/component/topology/container/dynamic/TriangleSetGeometryAlgorithms.h>
 
+
 namespace sofa::component
 {
 using sofa::component::topology::container::dynamic::TriangleSetTopologyModifier;
@@ -46,14 +47,14 @@ class TearingAlgorithms
 	using Edge = sofa::core::topology::BaseMeshTopology::Edge;
 	using Triangle = sofa::core::topology::BaseMeshTopology::Triangle;
 	using VecIds = sofa::type::vector<Index>;
+	
 
 	TearingAlgorithms(sofa::core::topology::BaseMeshTopology* _topology,
 		TriangleSetTopologyModifier* _modifier,
 		TriangleSetGeometryAlgorithms<DataTypes>* _triangleGeo);
 
 	virtual ~TearingAlgorithms();
 
-
 	/// <summary>
 	/// compute fracture path intersection point and cut through them
 	/// </summary>
@@ -158,6 +159,7 @@ class TearingAlgorithms
 
 	/// path created by algoFracturePath
 	sofa::type::vector<Coord> m_fracturePath;
+
 };
 
 

src/Tearing/TearingAlgorithms.inl

@@ -24,6 +24,7 @@
 #pragma once
 #include <Tearing/TearingAlgorithms.h>
 
+
 namespace sofa::component
 {
 

src/Tearing/TearingEngine.inl

@@ -96,14 +96,12 @@ inline bool TearingEngine<DataTypes>::computeIntersectionNeighborTriangle(Coord
 template<class DataTypes>
 inline bool TearingEngine<DataTypes>::computeEndPointsNeighboringTriangles(Coord Pa, Coord direction, 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;
 
@@ -151,12 +149,6 @@ void TearingEngine<DataTypes>::algoFracturePath()
     Coord Pb;
     Coord Pc;
 
-    if (this->d_fractureMaxLength.getValue())
-        this->computeEndPoints(Pa, principalStressDirection, Pb, Pc);
-    else if (!(computeEndPointsNeighboringTriangles(Pa, principalStressDirection, Pb, Pc)))
-        return;
-
-
     this->m_tearingAlgo->algoFracturePath(Pa, indexA, Pb, Pc, m_maxStressTriangleIndex, principalStressDirection, d_input_positions.getValue());
     m_maxStressTriangleIndex = InvalidID;
 
@@ -168,7 +160,7 @@ void TearingEngine<DataTypes>::algoFracturePath()
 template <class DataTypes>
 void TearingEngine<DataTypes>::computeFracturePath()
 {
-    if (!this->d_fractureMaxLength.getValue() && m_maxStressTriangleIndex != InvalidID)
+    if (m_maxStressTriangleIndex != InvalidID)
     {
         //Recording the endpoints of the fracture segment
         helper::ReadAccessor< Data<VecCoord> > x(d_input_positions);
@@ -179,11 +171,20 @@ void TearingEngine<DataTypes>::computeFracturePath()
         Coord Pb, Pc;
         fractureSegmentEndpoints.clear();
 
-        if (computeEndPointsNeighboringTriangles(Pa, principalStressDirection, Pb, Pc))
+        if (this->d_fractureMaxLength.getValue() == 0.0) {
+            computeEndPointsNeighboringTriangles(Pa, principalStressDirection, Pb, Pc);
+        }
+        else 
         {
-            fractureSegmentEndpoints.push_back(Pb);
-            fractureSegmentEndpoints.push_back(Pc);
+            this->computeEndPoints(Pa, principalStressDirection, Pb, Pc);
         }
+
+        fractureSegmentEndpoints.push_back(Pb);
+        fractureSegmentEndpoints.push_back(Pc);
+
+        this->m_stepCounter++;
+
+        //this->m_tearingAlgo->computeFracturePath(Pa, m_maxStressTriangleIndex, Pb, Pc);
     }
 }