Merge pull request #8 from ORNL/critical_stretch_condition_fix

Fixed condition for critical stretch criterion
ORNL · Aug 30, 2023 · c809608 · c809608
2 parents e58dd57 + 45aba40
commit c809608
Showing 1 changed file with 6 additions and 3 deletions.
diff --git a/Source/BondBreaking.m b/Source/BondBreaking.m
@@ -41,6 +41,9 @@
 
 function [u_NA] = BondBreaking(xx,yy,v,w,so,u_NA,r_hat_NA,x_hat_NA,y_hat_NA,mask_nofail)
 
+    % Tolerance
+    tol = 1E-15;
+
     % Number of nodes
     Nnodes = length(xx);
 
@@ -100,12 +103,12 @@
                 rk2 = rxk^2 + ryk^2;       % current bond length squared
 
                 % Note: the critical stretch condition is:
-                %       s := (sqrt(rk2) - Rk)/Rk > so
+                %       s := (sqrt(rk2) - Rk)/Rk >= so
                 %       This is equivalent to the condition
-                %             rk2 > Rk^2 * (so+1)^2
+                %             rk2 >= Rk^2 * (so+1)^2
                 %       because s = sqrt(rk2)/Rk - 1
 
-                if rk2  > ((so+1)*Rk)^2
+                if rk2  > ((so+1)*Rk)^2 - tol
 
                     % Remove cell uk from neighbor list of node ui
                     u_NA(ui,z) = 0;