%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof,emat,bc,F_ext,mat,ndim,nel,node,ndof,nip,nlod] = ex_unity; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol  = 0.00000001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material parameters for density growth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 1;
nue  = 0.3;
rho0 = 1.0;
psi0 = 1.0;
expm = 3.0;
expn = 2.0;
dt   = 1.0;

mat = [emod,nue,rho0,psi0,expm,expn,dt];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material parameters for volume growth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 1.0;
nue  = 0.3;
kt   = 2.0;
kc   = 2.0;
mt   = 2.0;
mc   = 2.0;
tt   = 2.0;
tc   = 0.5;
dt   = 1.0;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate node points q0 and connectivity edof
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xbox(1) = 0.0;
xbox(2) = 1.0;
ybox(1) = 0.0;
ybox(2) = 1.0;

nx =  2;
ny =  1;
nx =  2;
ny =  2;

[q0,edof] = mesh_sqr(xbox,ybox,nx,ny);  

[nel, sizen] = size(edof);
[ndof,sizen] = size(q0);
 for ie=1:nel emat(ie) = 1; end;
 node        = ndof/2;
 nip  = 4;
 ndim = 2;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dirichlet boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xmin = min(q0(1:2:ndof));
xmax = max(q0(1:2:ndof));
ymin = min(q0(2:2:ndof));
ymax = max(q0(2:2:ndof));

dy = ymax / 100;

nb = 1;

lmin=xmax/2-xmax/nx;
lmax=xmax/2+xmax/nx;

for in=2:2:ndof
  if((q0(in-1)<=xmin+tol)&(q0(in)<=ymin+tol))
    bc(nb,1) = in;
    bc(nb,2) = 0;
    nb = nb+1;
  end  
  if((q0(in-1)>=xmax-tol)&(q0(in)<=ymin+tol))
    bc(nb,1) = in;
    bc(nb,2) = 0;
    nb = nb+1;
  end  
  if(((q0(in-1)>=lmin-tol)&(q0(in-1)<=lmax+tol))&(q0(in)>=ymax-tol))
    bc(nb,1) = in;
    bc(nb,2) = -dy;
    nb = nb+1;
  end   
  if(((q0(in-1)>=xmax/2-tol)&(q0(in-1)<=xmax/2+tol))&(q0(in)<=ymin+tol))
    bc(nb,1) = in-1;
    bc(nb,2) = 0;
    nb = nb+1;
  end   
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% neumann boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_ext =  zeros(ndof,1);

%bc(1,1)= 1; bc(1,2)=0;
%bc(2,1)= 2; bc(2,2)=0;
%bc(3,1)= 3; bc(3,2)=0;
%bc(4,1)= 6; bc(4,2)=0;
%bc(5,1)=10; bc(5,2)=0;

%bc(6,1)= 9; bc(6,2)=0.0125;
%bc(7,1)=11; bc(7,2)=0.025;

bc(1,1)= 1; bc(1,2)=0;
bc(2,1)= 2; bc(2,2)=0;
bc(3,1)= 3; bc(3,2)=0;
bc(4,1)= 5; bc(4,2)=0;
bc(5,1)= 8; bc(5,2)=0;
bc(6,1)=14; bc(6,2)=0;

bc(7,1)=15; bc(7,2)=0.0125;
bc(8,1)=17; bc(8,2)=0.025;

%F_ext(13) = -0.002;

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