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

emod = 500;
nue  = 0.45;
rho0 = 1.00;
psi0 = 0.01;
expm = 3.0;
expn = 2.0;
dt   = 2.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

nx =  4;  ny =  4;  nz =   4;
lx =  1;  ly =  1;  lz =   1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% generate node points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ny = nx;

n = 0;
nd= 0;
for iz=0:nz
  for ix=0:nx
    for iy=0:ny
      nd = nd+1;
      n  = n +3;
      x0(nd,1)  = ix * lx/nx;
      y0(nd,1)  = iy * ly/ny;
      z0(nd,1)  = iz * lz/nz;
      q0(n-2,1) = ix * lx/nx;
      q0(n-1,1) = iy * ly/ny;
      q0(n  ,1) = iz * lz/nz;
    end
  end
end
[ndofc,sizen] = size(q0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% generate cubes (with additional center nodes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elem = 0;
for iz=1:nz
  for ix=1:nx
    for iy=1:ny
      elem = elem + 1;
      azu = 3*( iz-1)*(nx+1) * (nx+1);
      azo = 3*( iz  )*(nx+1) * (nx+1);
      nlu = 3*((ix-1)*(nx+1) + (iy-1))+1 +azu;
      nru = 3*( ix   *(nx+1) + (iy-1))+1 +azu;
      nlo = 3*((ix-1)*(nx+1) + (iy-1))+1 +azo;
      nro = 3*( ix   *(nx+1) + (iy-1))+1 +azo;
      
      ebox(elem,:) = [elem  nlu  nlu+1 nlu+2 nru nru+1 nru+2 nru+3 nru+4 nru+5 nlu+3 nlu+4 nlu+5 ...
                            nlo  nlo+1 nlo+2 nro nro+1 nro+2 nro+3 nro+4 nro+5 nlo+3 nlo+4 nlo+5];                                  
      n  = n +3;
      q0(n-2,1) = (q0(nlu  )+q0(nru  )+q0(nru+3)+q0(nlu+3)+q0(nlo  )+q0(nro  )+q0(nro+3)+q0(nlo+3))/8;
      q0(n-1,1) = (q0(nlu+1)+q0(nru+1)+q0(nru+4)+q0(nlu+4)+q0(nlo+1)+q0(nro+1)+q0(nro+4)+q0(nlo+4))/8;
      q0(n  ,1) = (q0(nlu+2)+q0(nru+2)+q0(nru+5)+q0(nlu+5)+q0(nlo+2)+q0(nro+2)+q0(nro+5)+q0(nlo+5))/8;
      
    end
  end
end  

index=[4;7;10;13;16;19;22;25];
enod=ebox(:,index)/3;
[nelc, sizen] = size(ebox);
[ndofp,sizen] = size(q0);

disp(['*** step [1] *** ',num2str(nelc), ' cubes with ', ...
                          num2str(ndofc/3),' nodes ***']);      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% generate tets (with additional dual nodes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ntet = 0;

for ie=1:elem
  ntet = ntet+12;
  center = nd +ie;
  etet(ntet-11,:) =[enod(ie,1) enod(ie,2) enod(ie,3) center];
  etet(ntet-10,:) =[enod(ie,1) enod(ie,3) enod(ie,4) center]; 
  etet(ntet-9 ,:) =[enod(ie,1) enod(ie,5) enod(ie,6) center];
  etet(ntet-8 ,:) =[enod(ie,1) enod(ie,6) enod(ie,2) center];
  etet(ntet-7 ,:) =[enod(ie,1) enod(ie,4) enod(ie,8) center];
  etet(ntet-6 ,:) =[enod(ie,1) enod(ie,8) enod(ie,5) center];
  etet(ntet-5 ,:) =[enod(ie,7) enod(ie,3) enod(ie,2) center];
  etet(ntet-4 ,:) =[enod(ie,7) enod(ie,2) enod(ie,6) center];
  etet(ntet-3 ,:) =[enod(ie,7) enod(ie,8) enod(ie,4) center];
  etet(ntet-2 ,:) =[enod(ie,7) enod(ie,4) enod(ie,3) center];
  etet(ntet-1 ,:) =[enod(ie,7) enod(ie,6) enod(ie,5) center];
  etet(ntet-0 ,:) =[enod(ie,7) enod(ie,5) enod(ie,8) center];
end

for it=1:ntet
  edof(it,:)=[it (etet(it,1)*3-2)  (etet(it,1)*3-1)  (etet(it,1)*3) ...
                 (etet(it,2)*3-2)  (etet(it,2)*3-1)  (etet(it,2)*3) ...
                 (etet(it,3)*3-2)  (etet(it,3)*3-1)  (etet(it,3)*3) ...
                 (etet(it,4)*3-2)  (etet(it,4)*3-1)  (etet(it,4)*3)];
end

[nel, sizen] = size(edof);
[ndof,sizen] = size(q0);
 for ie=1:nel emat(ie) = 1; end;
 
 node        = ndof/3;
 nip         = 1;
 ndim        = 3;

e_mat  = extr_dof(edof,q0);      
disp(['*** step [2] *** ',num2str(nel), ' primal tets with ', ...
                          num2str(ndofp/3),' nodes ***']);      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dirichlet boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

dz = zmax / 100;

nb = 1;

lmin= 1/4*xmax;
lmax= 3/4*xmax;

lmin=   0*xmax;
lmax= 2/4*xmax;

F_ext    =  zeros(ndof,1);

for in=3:3:ndof
  if(q0(in)<=zmin)    % unten 
%    bc(nb,1) = in-2;
%    bc(nb,2) = 0;
%    nb = nb+1;
%    bc(nb,1) = in-1;
%    bc(nb,2) = 0;
%    nb = nb+1;
    bc(nb,1) = in;
    bc(nb,2) = 0;
    nb = nb+1;
  end  

%  if(q0(in)>=zmax)    % oben 
%    bc(nb,1) = in-2;
%    bc(nb,2) = 0.01;
%    nb = nb+1;
%    bc(nb,1) = in-1;
%    bc(nb,2) = 0;
%    nb = nb+1;
%    bc(nb,1) = in;
%    bc(nb,2) = 0;
%    nb = nb+1;
%  end  

  if(q0(in-2)==xmin)  % links
    bc(nb,1) = in-2;
    bc(nb,2) = 0;
    nb = nb+1;
  end  

  if(q0(in-2)==xmax)  % rechts
    bc(nb,1) = in-2;
    bc(nb,2) = 0;
    nb = nb+1;
  end  

  if(q0(in-1)==ymin)  % hinten
    bc(nb,1) = in-1;
    bc(nb,2) = 0;
    nb = nb+1;
  end  
  
  if(q0(in-1)==ymax)  % vorne
    bc(nb,1) = in-1;
    bc(nb,2) = 0;
    nb = nb+1;
  end  

%  punch
  if(((q0(in-2)>=lmin)&(q0(in-2)<=lmax)) ...
    &((q0(in-1)>=lmin)&(q0(in-1)<=lmax)) ...
                      &(q0(in  )>=zmax))
    bc(nb,1) = in-2;
    bc(nb,2) = 0;
    nb = nb+1;
    bc(nb,1) = in-1;
    bc(nb,2) = 0;
    nb = nb+1;
    bc(nb,1) = in;
    bc(nb,2) = -dz;
    nb = nb+1;
%	F_ext(in,1)=-0.1;
  end   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% neumann boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%F_ext    =  zeros(ndof,1);
nlod     = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nmat,ngrw] = size(mat);
if (ngrw==7); disp('*** density growth ***'); end;  % density
if (ngrw==9); disp('*** volume growth ***' ); end;  % volume
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

