%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof,emat,bc,F_ext,mat,ndim,nel,node,ndof,nip,nlod] = ex_cylin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 500;
nue  = 0.3;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l    =   8.0;             % length
a    =   4.0;             % cross section
nel1 =     4;             % number of elements in plane - scaled
nel2 =     8;             % number of elements in longitudinal direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol  =  1e-8;
ndim = 3;
nip  = 8;

dp  =  1;
dp  =  0;
du  = l/20;
%du  =  0;
nls = 10;

blck = 1;
r    = sqrt(a/pi);

nip  = 8;
nep  = nel1*nel1*3*4;
nel  = nep*nel2;
node = ((nel1*2+1)*(nel1*2+1)+nel1*nel1*4*2)*(nel2+1);
ndof = ndim * node;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% coordinates %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q0 = zeros(ndim*node,1);
nn = 0;   

for iz = 0:nel2                  
%%% generation of inner block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    for ix = 0:2*nel1
	for iy = 0:2*nel1
	    x = -r/2.0+ix*r/2.0/nel1;
		y = -r/2.0+iy*r/2.0/nel1;
		z = l/nel2*iz; 
	    nn = nn +ndim;
        q0(nn-2,1) = x;
        q0(nn-1,1) = y;
        q0(nn  ,1) = z;
    end
	end  
      
%%% generation of outer block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    for ix = 1:nel1
    for iy = 1:(8*nel1)
		alpha =  pi*(iy-1)/4.0/nel1 - pi/4.0;
		if  cos(alpha) == 0
		    b = abs(r/2.0/sin(alpha));
		else 
			if sin(alpha) == 0
			   b = abs(r/2.0/cos(alpha));
			else		     	
			   b = min(abs(r/2.0/cos(alpha)),abs(r/2.0/sin(alpha)));
			end	
		end
		dl=  b+ix*(r-b)/nel1;
        x = -cos(alpha)*dl;
		y =  sin(alpha)*dl;
		z = l/nel2*iz;
	    nn = nn +ndim;
        q0(nn-2,1) = x;
        q0(nn-1,1) = y;
        q0(nn  ,1) = z;
    end
	end   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% connectivity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
enod = zeros(nel,8);
edof = zeros(nel,8*ndim+1);
ee   = 0;

for iz=1:nel2                   
	nodl=((nel1*2+1)*(nel1*2+1)+nel1*nel1*4*2)*(iz-1)+1;
	nodu=((nel1*2+1)*(nel1*2+1)+nel1*nel1*4*2)*(iz)  +1;         
	dnod= (nel1*2+1)*(nel1*2+1);
         
! ****** generation of inner block *********************************** !     
    for ix=1:(2*nel1)
	for iy=1:(2*nel1)
	    ee = ee + 1;
		n1 = nodu+(2*nel1+1)*(ix-1)+iy-1;
		n2 = nodu+(2*nel1+1)*(ix  )+iy-1;
		n3 = nodu+(2*nel1+1)*(ix  )+iy;
		n4 = nodu+(2*nel1+1)*(ix-1)+iy;
		n5 = nodl+(2*nel1+1)*(ix-1)+iy-1;
		n6 = nodl+(2*nel1+1)*(ix  )+iy-1;
		n7 = nodl+(2*nel1+1)*(ix  )+iy;
		n8 = nodl+(2*nel1+1)*(ix-1)+iy;
		enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
    end
	end
         
! ****** generation of connecting block ****************************** !     
    for ix=1:nel1*2 
	    ee = ee + 1;
        n1 = nodu-1+ix;
		n2 = nodu-1+ix     +1;
		n3 = nodu-1+ix+dnod+1;
		n4 = nodu-1+ix+dnod;
		n5 = nodl-1+ix;
		n6 = nodl-1+ix     +1;
		n7 = nodl-1+ix+dnod+1;
		n8 = nodl-1+ix+dnod;
		enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
    end

	for ix=1:nel1*2 
	    ee = ee + 1;
        n1 = nodu-1+ix*(nel1*2+1);
        n2 = nodu-1+(ix+1)*(nel1*2+1);
        n3 = nodu-1+nel1*2+ix+dnod+1;
        n4 = nodu-1+nel1*2+ix+dnod;
        n5 = nodl-1+ix*(nel1*2+1);
        n6 = nodl-1+(ix+1)*(nel1*2+1);
        n7 = nodl-1+nel1*2+ix+dnod+1;
        n8 = nodl-1+nel1*2+ix+dnod;
		enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
    end

    for ix=1:nel1*2 
	    ee = ee + 1;
        n1 = nodu+dnod-ix;
        n2 = nodu+dnod-ix-1;
        n3 = nodu-1+nel1*4+ix+dnod+1;
        n4 = nodu-1+nel1*4+ix+dnod;
        n5 = nodl+dnod-ix;
        n6 = nodl+dnod-ix-1;
        n7 = nodl-1+nel1*4+ix+dnod+1;
        n8 = nodl-1+nel1*4+ix+dnod;
		enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
    end

         for ix=1:nel1*2-1 
	     ee = ee + 1;
         n1 = nodu+dnod-ix*(nel1*2+1);
         n2 = nodu+dnod-(ix+1)*(nel1*2+1);
         n3 = nodu-1+nel1*6+ix+dnod+1;
         n4 = nodu-1+nel1*6+ix+dnod;
         n5 = nodl+dnod-ix*(nel1*2+1);
         n6 = nodl+dnod-(ix+1)*(nel1*2+1);
         n7 = nodl-1+nel1*6+ix+dnod+1;
         n8 = nodl-1+nel1*6+ix+dnod;
		 enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
    end
 	     ee = ee + 1;
         n1 = nodu+dnod-(nel1*2)*(nel1*2+1);
         n2 = nodu+dnod-(nel1*2+1)*(nel1*2+1);
         n3 = nodu+dnod;
         n4 = nodu-1+nel1*6+(nel1*2)+dnod;
         n5 = nodl+dnod-(nel1*2)*(nel1*2+1);
         n6 = nodl+dnod-(nel1*2+1)*(nel1*2+1);
         n7 = nodl+dnod;
         n8 = nodl-1+nel1*6+(nel1*2)+dnod;
 		 enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
         
            
! ****** generation of outer block *********************************** !     
         for ix=1:nel1-1
         for iy=1:(8*nel1-1)
 	     ee = ee + 1;
         n1 = nodu+dnod+(ix-1)*(4*2*nel1)+iy-1;
         n2 = nodu+dnod+(ix-1)*(4*2*nel1)+iy;
         n3 = nodu+dnod+ ix   *(4*2*nel1)+iy;
         n4 = nodu+dnod+ ix   *(4*2*nel1)+iy-1;
         n5 = nodl+dnod+(ix-1)*(4*2*nel1)+iy-1;
         n6 = nodl+dnod+(ix-1)*(4*2*nel1)+iy;
         n7 = nodl+dnod+ ix   *(4*2*nel1)+iy;
         n8 = nodl+dnod+ ix   *(4*2*nel1)+iy-1;
 		 enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
   end
 	     ee = ee + 1;
		 n1 = nodu+dnod+(ix-1)*(4*2*nel1)+4*2*nel1-1;
         n2 = nodu+dnod+(ix-1)*(4*2*nel1);
         n3 = nodu+dnod+(ix-1)*(4*2*nel1)+4*2*nel1;
         n4 = nodu+dnod+ ix   *(4*2*nel1)+4*2*nel1-1;
         n5 = nodl+dnod+(ix-1)*(4*2*nel1)+4*2*nel1-1;
         n6 = nodl+dnod+(ix-1)*(4*2*nel1);
         n7 = nodl+dnod+(ix-1)*(4*2*nel1)+4*2*nel1;
         n8 = nodl+dnod+ ix   *(4*2*nel1)+4*2*nel1-1;             
 		 enod(ee,:) = [n1  n2  n3  n4  n5  n6  n7  n8];
   end   
end

for ie = 1:nel
    edof(ie,:) = [ie, ndim*enod(ie,5)-2  ndim*enod(ie,5)-1  ndim*enod(ie,5) ...
					  ndim*enod(ie,6)-2  ndim*enod(ie,6)-1  ndim*enod(ie,6) ...
					  ndim*enod(ie,7)-2  ndim*enod(ie,7)-1  ndim*enod(ie,7) ...
					  ndim*enod(ie,8)-2  ndim*enod(ie,8)-1  ndim*enod(ie,8) ...
					  ndim*enod(ie,1)-2  ndim*enod(ie,1)-1  ndim*enod(ie,1) ...
                      ndim*enod(ie,2)-2  ndim*enod(ie,2)-1  ndim*enod(ie,2) ...
	  				  ndim*enod(ie,3)-2  ndim*enod(ie,3)-1  ndim*enod(ie,3) ...
					  ndim*enod(ie,4)-2  ndim*enod(ie,4)-1  ndim*enod(ie,4) ];
end
for ie=1:nel emat(ie) = 1; end;

%%% boundary conditions & loading %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% all boudary fixed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nb  = 0;
npl = (nel1*2+1)*(nel1*2+1)+nel1*nel1*4*2;

for ib = 1:npl
  nb = nb+1;
%  bc(nb,:) = [ndim*ib-2  0];
  nb = nb+1;
%  bc(nb,:) = [ndim*ib-1  0];
  nb = nb+1;
%  bc(nb,:) = [ndim*ib-0  0];
end
for ib = 1:npl
  nb = nb+1;
%  bc(nb,:) = [ndim*(node-npl)+ndim*ib-2  0];
  nb = nb+1;
%  bc(nb,:) = [ndim*(node-npl)+ndim*ib-1  0];
  nb = nb+1;
%  bc(nb,:) = [ndim*(node-npl)+ndim*ib-0  du];
end

%%% homogeneous problem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nb = 0;
for ib = 1:npl
  if(abs(q0(ndim*ib-2)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*ib-2  0];
  end
  if(abs(q0(ndim*ib-1)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*ib-1  0];
  end	
    nb = nb+1;
    bc(nb,:) = [ndim*ib-0  0];
end

for ib = 1:npl
  if(abs(q0(ndim*(node-npl)+ndim*ib-2)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*(node-npl)+ndim*ib-2  0];
  end
  if(abs(q0(ndim*(node-npl)+ndim*ib-1)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*(node-npl)+ndim*ib-1  0];
  end
    nb = nb+1;
    bc(nb,:) = [ndim*(node-npl)+ndim*ib-0  du];
end

%%% point loads %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nlod  = 5;
load  = 0.0;
F_ext = zeros(ndof,1); 
F_ext( 3)=load;
F_ext( 6)=load;
F_ext( 9)=load;
F_ext(12)=load;
F_ext = zeros(ndof,1); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nmat,ngrw] = size(mat);
if (ngrw==7); disp('*** density growth ***'); end;  % density
if (ngrw==9); disp('*** volume growth ***' ); end;  % volume
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
