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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material parameters for volume growth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 100;
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) = 8.0;
ybox(1) = 0.0;
ybox(2) = 1.0;

nx = 32;
ny =  4;

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

[nel, sizen] = size(edof);
[ndof,sizen] = size(q0);
 node        = ndof/2;
 for ie=1:nel emat(ie) = 1; end;
 
 nip  = 4;
 ndim = 2;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dirichlet boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fdof1 = 2*(ny+1)*(0     ) +2;
    fdof2 = 2*(ny+1)*(1     ) +2;
    fdof3 = 2*(ny+1)*(nx-1  ) +2;
    fdof4 = 2*(ny+1)*(nx    ) +2;
	fdof5 = 2*(ny+1)*(nx/2+0) +1;

    fdof0 = 2*(ny+1)*(0     ) +1;
    fdof9 = 2*(ny+1)*(nx    ) +1;

	fdof6 = 2*(ny+1)*(nx/2  ); 
	fdof7 = 2*(ny+1)*(nx/2+1);
	fdof8 = 2*(ny+1)*(nx/2+2); 
	
%    bc(1,1) = fdof1; bc(1,2) = 0;
%    bc(2,1) = fdof4; bc(2,2) = 0;   
%	bc(3,1) = fdof5; bc(3,2) = 0;
%	bc(4,1) = fdof7; bc(4,2) = -ybox(2)/50;
% 	bc(5,1) = fdof6; bc(5,2) = -ybox(2)/50;
% 	bc(6,1) = fdof8; bc(6,2) = -ybox(2)/50;
%%    bc(7,1) = fdof2; bc(7,2) = 0;   
%%	 bc(8,1) = fdof3; bc(8,2) = 0;
%%    bc(9,1) = fdof0; bc(9,2) = 0;
%%    bc(10,1) =fdof9; bc(10,2)= 0; 

    bc(1,1) =   6; bc(1,2) = 0;
	bc(2,1) = 159; bc(2,2) = 0;  
	bc(3,1) = 326; bc(3,2) = 0;  
	
%	bc(3,1) = 50; bc(3,2) = -ybox(2)/50;
%	bc(4,1) = 86; bc(4,2) = 0;  

			
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% neumann boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fp = -0.002;

F_ext =  zeros(ndof,1);
F_ext(10) = Fp/2;
F_ext(20) = Fp;
F_ext(30) = Fp;
F_ext(40) = Fp;
F_ext(50) = Fp;
F_ext(60) = Fp;
F_ext(70) = Fp;
F_ext(80) = Fp;
F_ext(90) = Fp;
F_ext(100) = Fp;
F_ext(110) = Fp;
F_ext(120) = Fp;
F_ext(130) = Fp;
F_ext(140) = Fp;
F_ext(150) = Fp;
F_ext(160) = Fp;
F_ext(170) = Fp;
F_ext(180) = Fp;
F_ext(190) = Fp;
F_ext(200) = Fp;
F_ext(210) = Fp;
F_ext(220) = Fp;
F_ext(230) = Fp;
F_ext(240) = Fp;
F_ext(250) = Fp;
F_ext(260) = Fp;
F_ext(270) = Fp;
F_ext(280) = Fp;
F_ext(290) = Fp;
F_ext(300) = Fp;
F_ext(310) = Fp;
F_ext(320) = Fp;
F_ext(330) = Fp/2;

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