%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof,emat,bc,F_ext,mat,ndim,nel,node,ndof,nip,nlod] = ex_femur
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 500;
nue  = 0.2;
rho0 = 1.00;
psi0 = 0.01;
expm = 3.0;
expn = 2.0;
dt   = 5.0;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material parameters for volume growth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 500;
nue  = 0.2;
kt   = 1.0;
kc   = 1.0;
mt   = 1.2;
mc   = 1.2;
tt   = 1.1;
tc   = 0.9;
dt   = 0.001;

% mat = [emod,nue,kt,kc,mt,mc,tt,tc,dt];

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

[nel,sizen]  = size(edof);
[ndof,sizen] = size(q0);
 
 for ie=1:nel emat(ie) = 1; end;
 
 node = ndof/2;
 nip  = 4;
 ndim = 2;
 tol  =  1e-8;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dirichlet boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
  bc(  1,1) = 1;     bc(1,2) = 0;    % fix node    1 in x-dir
  bc(  2,1) = 2;     bc(2,2) = 0;    % fix node    1 in y-dir
for i=2:15
  bc(i+1,1) = 2*i;   bc(i+1,2) = 0;  % fix node 2-15 in x-dir
end
			
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% neumann boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_ext = zeros(ndof,1);
F_ext(372*2-1)= 0.5496;  F_ext(372*2)= 1.3517; % load case 1,2,3
F_ext(692*2-1)= 0.2997;  F_ext(692*2)=-1.1185; % load case 2
F_ext(720*2-1)=-1.2792;  F_ext(720*2)=-0.8628; % load case 3
F_ext(723*2-1)=-0.9424;  F_ext(723*2)=-2.1167; % load case 1

nlod = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nmat,ngrw] = size(mat);
if (ngrw==7); disp('*** density growth ***'); end;  % density
if (ngrw==9); disp('*** volume growth ***' ); end;  % volume
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%