%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof,emat,bc,F_ext,mat,ndim,nel,node,ndof,nip,nlod] = ex_tube2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% material parameters for volume growth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 3.0;
nue  = 0.3;
kt   = 0.5;
kc   = 0.25;
mt   = 4.0;
mc   = 5.0;
tt   = 1.5;
tc   = 0.5;
dt   = 0.5;

mat = [emod,nue,kt,kc,mt,mc,tt,tc,dt];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l   =    2.0;      % length      10.0
ra  =    1.0;      % inner radius 1.3
ri  =    0.5;      % outer radius 1.0
nez =      8;      % number of elements in longitudinal direction %21
ner =      4;      % number of elements in radial direction
nep =     16;      % number of elements in circumferential direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol  = 1e-8;
ndim = 3;
nip  = 8;
nel = nez * ner * nep;
node= (nez+1)*(ner+1)*nep;
ndof = ndim*node;
%%% coordinates %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q0 = zeros(ndim*node,1);
nn = 0;   

delta_z =       l / nez;
delta_r = (ra-ri) / ner;
delta_t =    2*pi / nep;

for iz = 0:nez
  z = iz * delta_z;
  for ir = 0:ner
    r = ri + ir * delta_r;
	for ip = 0:(nep-1)
	  p = ip * delta_t; 
      nn = nn +ndim;
      q0(nn-2,1) = r*cos(p);
      q0(nn-1,1) = r*sin(p);
      q0(nn  ,1) = z;
	end
  end	  
end

%%% connectivity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
enod = zeros(nel,8);
edof = zeros(nel,8*ndim+1);
ee   = 0;
n1   = 0;
 
for iz = 1:nez
  for ir = 1:ner
	for ip = (1:nep-1)
	  ee = ee + 1; 
      n1 = n1 + 1;
	  n2 = n1 + nep*(ner+1);
	  n3 = n1 + nep;
	  n4 = n3 + nep*(ner+1);
	  enod(ee,:) = [n1  n1+1  n2+1  n2  n3  n3+1  n4+1  n4];
	end
	ee = ee+1;
    n1 = n1 + 1;
	n2 = n1 + nep*(ner+1);
	n3 = n1 + nep;
    n4 = n3 + nep*(ner+1);
    enod(ee,:) = [n1  n1-nep+1  n2-nep+1  n2  n3  n3-nep+1  n4-nep+1  n4];
  end	  
  n1 = n1+nep;
end

for ie = 1:nel
  edof(ie,:) = [ie, 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) ...
					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)];
end
for ie=1:nel emat(ie) = 1; end;
 
%%% boundary conditions & loading %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% homogeneous problem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
du = -l/10;
nb = 0;

for ib = 1:(nep*(ner+1))
  if(abs(q0(ndim*(node-nep*(ner+1))+ndim*ib-2)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*ib-2  0];
  end
  if(abs(q0(ndim*(node-nep*(ner+1))+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:(nep*(ner+1))
  if(abs(q0(ndim*(node-nep*(ner+1))+ndim*ib-2)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*(node-nep*(ner+1))+ndim*ib-2  0];
  end
  if(abs(q0(ndim*(node-nep*(ner+1))+ndim*ib-1)-0.0)<tol)
    nb = nb+1;
    bc(nb,:) = [ndim*(node-nep*(ner+1))+ndim*ib-1  0];
  end
  nb = nb+1;
  bc(nb,:) = [ndim*(node-nep*(ner+1))+ndim*ib-0  du];
end

%%% point loads %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load  = 0.0;
F_ext = zeros(ndof,1);  
nlod  = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['*** generated ',num2str(nel),' *** elems with ',num2str(ndof),' dofs ***']);      
[nmat,ngrw] = size(mat);
if (ngrw==7); disp('*** density growth ***'); end;  % density
if (ngrw==9); disp('*** volume growth ***' ); end;  % volume
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


