# Aleksandar Donev, December 2002
# Models for jamming in hard-sphere packings
# Testing for jamming in sphere packings:

fix {id in IDs} dD[id]; # No attempt to grow the particles
let {id in IDs} D_perturbed[id]:=D_relaxation[id]*D_perturbed[id]; # Try some relaxation	

# We must adjust the diameters so there is no overlap:

if 1 then {
   # We must ensure that there is no overlap
	for {id_1 in 1..n_ids} for {id_2 in id_1..n_ids}
	{
	   let overlap[id_1, id_2]:=min{a in ARCS: 
	      min(ID[H[a]],ID[T[a]])=id_1 and max(ID[H[a]],ID[T[a]])=id_2} l[a]/l_strut[a];
	   let overlap[id_2, id_1]:=overlap[id_1, id_2];
      print "Overlap: ", id_1, id_2, overlap[id_1, id_2];
	}
	let {id in IDs} D_relaxation[id]:=min{id_1 in IDs, id_2 in IDs} overlap[id_1,id_2];
	print "Relaxing diameters by:", {id in IDs} 1-D_relaxation[id];
	let {id in IDs} D_perturbed[id]:=D_relaxation[id]*D_perturbed[id];
	let {a in ARCS} l_strut[a]:=(D_perturbed[ID[H[a]]]+D_perturbed[ID[T[a]]])/2;
}	
	
# In troublesome cases:
print "BAD ARCS: ", {a in ARCS: (l_strut[a]-l[a])/l_strut[a]>1E-6} (a,(l_strut[a]-l[a])/l_strut[a]);

if n_frames=1 then let n_frames:=2*n_unjamming_loads*n_jamming_LPs;

objective Work;
print "Testing with relative jamming tolerance: ", displacement_tolerance;

for {load in 1..n_unjamming_loads} 
{
let {i in NODES, k in DIMS} b[i,k]:=Uniform(-1,1);
   
for {iteration in 1..2} 
{
   if iteration=2 then let {i in NODES, k in DIMS} b[i,k]:=-b[i,k];
      # Reverse the loading direction

   # Solve several LP's iteratively to make sure we explore nonlinearities!
   for {jamming_LP in 1..n_jamming_LPs}
   {      
	   let n_displaced:=0;
      let unjammed:=0;

      # Remove rattling clusters (for this load):	   	   
	   for {subpacking_attempt in 1..n_subpacking_attempts} 
	   {	      
	      print "Now at Load #: ", load, " of sign: ", iteration, 
	         " Unjamming attempt #: ", subpacking_attempt,
	         " Jamming LP #: ", jamming_LP;
	      solve; # Solve the LP
		   print "Average and maximal relative displacement of spheres: ",
			  (sum {i in NODES} (sqrt( sum{k in DIMS} (dr_scaling*dr[i,k])^2))/D_perturbed[ID[i]])/n_active,
			  (max {i in NODES} (sqrt( sum{k in DIMS} (dr_scaling*dr[i,k])^2))/D_perturbed[ID[i]]);
	      
			if test_subpacking then 
			{  
			   # We now remove the rattling spheres:
	         let n_rattling:=card({i in NODES: 
	            sqrt( sum{k in DIMS} dr[i,k]^2 ) > displacement_tolerance*D_perturbed[ID[i]]});	
				if(n_rattling=0) then {
				   print "Load seems supported by subpacking with ", n_active, " spheres!";
				   break; # This load is supported by subpacking
				}
	            		   
				let {i in NODES: sqrt( sum{k in DIMS} dr[i,k]^2) >=
				   displacement_tolerance*D_perturbed[ID[i]]} active_nodes[i]:=0; 
			   include Remove.mod; # Update other values as well		       	  
				print "Number of unjammed, jammed spheres:", n-n_active, n_active;				
				if (n_active=0 or m_active=0) then 
				{
				   print "There is no jammed subpacking!";
				   let unjammed:=1;
				   break;
				}			 		   
	      }
	   } # Finish iterative subpacking LP

      if unjammed then break; # No need to test non-linearities
      
		# Manually eliminate translations:
		let {k in DIMS} dr_average[k]:=(sum{i in NODES: active_nodes[i]} dr[i,k])/n_active;
		let {k in DIMS, i in NODES: active_nodes[i]} dr[i,k]:=dr[i,k]-dr_average[k];	         
	   include ScaleDR.mod;
	   include Displace.mod; # Update other dependent variables	            			   		

		print "Scaled strain is:";
		print {i in DIMS}: {j in DIMS} dr_scaling*strain[strain_indices[i,j]];

		let average_displacement:=
		   (sum {i in NODES} (sqrt( sum{k in DIMS} (r_perturbed[i,k]-r[i,k])^2))/D_perturbed[ID[i]])/n_active;
		let maximal_displacement:=
         (max {i in NODES} (sqrt( sum{k in DIMS} (r_perturbed[i,k]-r[i,k])^2))/D_perturbed[ID[i]]);		         
	   print "Average: ", average_displacement, 
	      " and maximal: ", maximal_displacement, " relative total displacement";
			
	   let n_displaced:=n_displaced+card({i in NODES: active_nodes[i] and
	      (sqrt( sum{k in DIMS} (r_perturbed[i,k]-r[i,k])^2 ) > displacement_tolerance*D_perturbed[ID[i]])});
	   print "Number of spheres displaced by load: ", n_displaced;
		if(n_displaced>0) then 
		{
		   let unjammed:=1;
		   break;
	   }	     
	   
	} # Finish iterative LP  
	if unjammed>0 then break;
} # Both +b and -b have been tried
if unjammed>0 then break;
} # All loads have been tried

if no_strain then
   {let plot_file := "Scenes/"&problem_name&".collective.unjamming.wrl";}
else
   {let plot_file := "Scenes/"&problem_name&".strict.unjamming.wrl";}       
include Plot.mod;   
   
let output_file := "Output/"&problem_name&".dr.dat";
if unjammed>0 then {
   print "SPHERE PACKING IS NOT JAMMED!!!";   
	   
   print "Writing unjamming motion to ", output_file;
	printf "%lf %d %d %d\n", l_min, n, 0, d > (output_file);
	print {i in DIMS, j in DIMS} L_perturbed[i,j] ,
	      {i in DIMS, j in DIMS} (sum {k in DIMS} dr_scaling*strain[strain_indices[i,k]]*L_perturbed[k,j]) 
	        > (output_file); # Output L and dL
	print {k in DIMS} 0.0 > (output_file);
	print {i in NODES}: if active_nodes[i] then "T" else "F",
	  {k in DIMS} r_perturbed[i,k], {k in DIMS} dr_scaling*dr[i,k] > (output_file);
	  # And also r and dr
   	  
}
else
{
   print "SPHERE PACKING IS JAMMED WITH ", n_active, " JAMMED SPHERES!!!";
}

# EOF