pro setup_slab
  fin='glass_10x10x10'
  readnew,fin,h,'HEAD'
  readnew,fin,xin,'POS'
  xin[*,*] /= h.BOXSIZE
  Nin = TOTAL(h.npart,/int)
  
  Ndup = 10L
  fout='glass_100x100x10'
  
  Np = Ndup^2 * Nin
  x=fltarr(3,Np)

  i=0L
  FOR ix=0,Ndup-1 DO BEGIN
     FOR iy=0,Ndup-1 DO BEGIN
        FOR j=0L,Nin-1 DO BEGIN
           x[0,i] = xin[0,j] + ix 
           x[1,i] = xin[1,j] + iy
           x[2,i] = xin[2,j]
           i++
        END
     END
  END
  
  v = fltarr(3,Np)
  u=fltarr(Np)
  m=fltarr(Np)
  id=lindgen(Np)

  h.BOXSIZE=1
  h.MASSARR[*]=0

  gamma=5./3.
  rho = 1.0
  P = 1.0
  
  u[*] = P/(gamma-1)/rho
  m[*] = rho *  (h.BOXSIZE^3 / Ndup) / Np

  h.npart[*]=0
  h.npart[0]=Np

  h.PARTTOTAL[*]=0
  h.PARTTOTAL[0]=Np

  dL = 1.0
  dv = 1.0
;Left moving Plate
   jj=WHERE(x[0,*] LT dL)
   id[jj] = 0
   v[1,jj] = -1 * dv
   u[jj] *= 1.05

;Right moving Plate
   jj=WHERE(x[0,*] GT Ndup-dL)
   id[jj] = 0
   v[1,jj] = +1 * dv
   u[jj] *= 1.05

; Bring positions to [0,1]
   x /= ndup
   
   write_head,fout,h
   add_block,fout,x,'POS '
   add_block,fout,v,'VEL '
   add_block,fout,id,'ID  '
   add_block,fout,m,'MASS'
   add_block,fout,u,'U   '

   jj=where(abs(x[2,*]-0.5) LT 0.1)
   plot,x[0,*],x[1,*],psym=3
   jj=where(id[*] EQ 0)
   oplot,x[0,jj],x[1,jj],psym=3,col=128
   
end