%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof,emat,bc,F_ext,mat,ndim,nel,node,ndof,nip,nlod] = ex_humr2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol  =  1e-8;

emod = 500;
nue  = 0.3;
rho0 = 1.00;
psi0 = 0.01;
psi0 = 1.00;
expm = 3.0;
expn = 2.0;
dt   = 1.0;
dt   = 10;

mat = [emod,nue,rho0,psi0,expm,expn,dt];
%material parameters
%emod  nue  psi0*  rho0* expm  expn   dt
%500  0.20  0.01   1.20  3.00  2.00  2.0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate node points q0 and connectivity edof
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[q0,edof]=in_humr2; 

[nel,sizen]  = size(edof);
[ndof,sizen] = size(q0);
 for ie=1:nel emat(ie) = 1; end;

 node = ndof/3;
 nip  = 1;
 ndim = 3;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dirichlet and neumann boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	   
%fixnod=[583];          % small elbow bone   
%fixnod=[ 91; 288; 633];% large elbow bone
%fixnod=[260; 610; 623];% shoulder
 fixnod=[ 91; 260; 288; 583; 610; 623; 633];

[nfix,one]=size(fixnod);		 		
nb = 1;
for i=1:nfix
    nd=fixnod(i);
    bc(nb,1) = nd*3-2;
    bc(nb,2) = 0;
    nb = nb+1;
    bc(nb,1) = nd*3-1;
    bc(nb,2) = 0;
    nb = nb+1;
    bc(nb,1) = nd*3;
    bc(nb,2) = 0;
    nb = nb+1;	
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
% phase III - ball impact - extended
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
F_ext = zeros(ndof,1);
nlod = 5;
F0 = 0.0005;

F_del = 1.600 * F0;
F_pec = 1.170 * F0;
F_lat = 1.380 * F0;
F_tri = 0.580 * F0;

% del
%F_ext(414*3-2,1) =+0.36 * F_del; 
%F_ext(414*3-1,1) =+0.89 * F_del;
%F_ext(414*3-0,1) =-0.27 * F_del;
% mpre aligned now
F_ext(414*3-2,1) =+0.26 * F_del;
F_ext(414*3-1,1) =+0.94 * F_del;
F_ext(414*3-0,1) =-0.22 * F_del;

% pec
F_ext(423*3-2,1) =+0.41 * F_pec;
F_ext(423*3-1,1) =+0.82 * F_pec;
F_ext(423*3-0,1) =-0.41 * F_pec;

% lat
F_ext(475*3-2,1) =-0.57 * F_lat;
F_ext(475*3-1,1) =+0.79 * F_lat;
F_ext(475*3-0,1) =-0.23 * F_lat;

% tri
%F_ext(240*3-2,1) =-0.36 * F_tri;  
%F_ext(240*3-1,1) =-0.89 * F_tri;
%F_ext(240*3-0,1) =+0.27 * F_tri;
% mpre aligned now
F_ext(240*3-2,1) =-0.26 * F_tri;
F_ext(240*3-1,1) =-0.94 * F_tri;
F_ext(240*3-0,1) =+0.22 * F_tri;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
% phase II - maximum external shoulder rotation - flexed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
F_ext = zeros(ndof,1);
nlod = 10;
F0 = 0.0005;
F0 = 0.1;

F_del = 1.410 * F0;
F_pec = 1.490 * F0;
F_lat = 1.275 * F0;
F_tri = 1.250 * F0;

% del
F_ext(414*3-2,1) =+0.46 * F_del;
F_ext(414*3-1,1) =+0.76 * F_del;
F_ext(414*3-0,1) =-0.46 * F_del;

% pec
F_ext(423*3-2,1) =+0.47 * F_pec;
F_ext(423*3-1,1) =+0.75 * F_pec;
F_ext(423*3-0,1) =-0.47 * F_pec;

% lat
F_ext(475*3-2,1) =+0.49 * F_lat;
F_ext(475*3-1,1) =+0.39 * F_lat;
F_ext(475*3-0,1) =-0.78 * F_lat;

% tri
F_ext(240*3-2,1) =-0.44 * F_tri;
F_ext(240*3-1,1) =-0.79 * F_tri;
F_ext(240*3-0,1) =+0.44 * F_tri;

% maybe turn off lat to bring in more rotation
% converges (too) fast in 5 steps loading and 15 steps iteration of dt = 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
% points to fix on Coarse Model:  
% 184,13,3,20,55,119,151,167,171,145        <-- or any 6 of these
%
% points to fix on Fine Model: 
% 180, 515, 20, 516, 55, 106, 167, 177, 514 <-- or any 6 of these
%
% points to apply force-- apply approx 200lbs at each node
%
% towards elbow    coarse: node #306   200lb * [0.1881  0.9407  -0.2822]
% towards shoulder coarse: node #310   200lb * [0.1826  0.9129  -0.3651]
%
% towards elbow    fine:   node #446   200lb * [0.1881  0.9407  -0.2822]
% towards shoulder fine:   node #453   200lb * [0.1826  0.9129  -0.3651]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
