# Aleksandar Donev, December 2002
# Models for jamming in hard-sphere packings
# Running displacement-formulation problems

reset;

model; include Packing.mod; include DR.mod;

# These are substituted by the shell for some automated runs:
let problem_name:="_problem_name";

let input_file:="Data/"&problem_name&".r.dat";
print "Reading file: ", input_file;
data (input_file);
if _read_extra then data ("Data/"&problem_name&".extra.dat"); # Read in additional information

let delta:=_delta; # Gap tolerance (important!)
let displacement_tolerance:=_displacement_tolerance;
let no_strain:=_no_strain; # Test for collective (1) or strict (0) jamming
include Packing.param; # Additional parameters

# option times 1;
option randseed _random_seed; # Generated by bash
option log_file ("Output/"&problem_name&".log");
option solver_msg 0;

if _use_kestrel then 
{
	option solver kestrel;
	option kestrel_options "solver=LOQO";
}
else
{
	option solver loqo;
}
option knitro_options "iprint=1 maxit=100"; # feasible=1 direct=1
option loqo_options ("iterlim=100 convex verbose="&solver_verbosity); # honor_bnds

if dilution>0 then {
	let {i in NODES} active_nodes[i]:= if Uniform(0,1)<dilution then 0 else 1;
	print "Randomly removed ", 
	      card({i in NODES: active_nodes[i]=0})," spheres!";
}

fix {k in DIMS} dr[1,k]; # Is it better constrain the total momentum to 0?

fix {k in DIMS, i in NODES: not active_nodes[i]} dr[i,k];
   # Freeze the first node and any other inactive nodes

if no_strain then {   
   drop strain_trace;
   fix strain;
}

display n, m+m_, (n-1)*d-d;
display n_active, m_active;
print "The average, max and min relative gaps are: ",
      (sum {a in ARCS} (l[a]-l_strut[a])/l_strut[a])/m_active, 
      (max {a in ARCS} (l[a]-l_strut[a])/l_strut[a]), 
      (min {a in ARCS} (l[a]-l_strut[a])/l_strut[a]);

# Test for isostaticity of contact network: 
print "The number of perfect (to within 1E-3) contacts is ", 
   card({a in ARCS: (l[a]-l_strut[a])/l_strut[a]<1E-3});
                    
# EOF