program cube
*             =====
      parameter (LAYOUT=51)
* ... Fortran statement functions,
* ... defining the shape functions of
* ... a Finite Element 8-node cube:
      s1(x,y,z)=(0.5-x)*(0.5-y)*(0.5-z)
      s2(x,y,z)=(0.5+x)*(0.5-y)*(0.5-z)
      s3(x,y,z)=(0.5-x)*(0.5+y)*(0.5-z)
      s4(x,y,z)=(0.5+x)*(0.5+y)*(0.5-z)
      s5(x,y,z)=(0.5-x)*(0.5-y)*(0.5+z)
      s6(x,y,z)=(0.5+x)*(0.5-y)*(0.5+z)
      s7(x,y,z)=(0.5-x)*(0.5+y)*(0.5+z)
      s8(x,y,z)=(0.5+x)*(0.5+y)*(0.5+z)
*
* ... 51 x 51 x 51 grid:
      do 10 i=1,LAYOUT
      xp=5.*(float(i-1)/float(LAYOUT-1)-0.5)
      do 15 j=1,LAYOUT
      yp=5.*(float(j-1)/float(LAYOUT-1)-0.5)
      do 20 k=1,LAYOUT
      zp=5.*(float(k-1)/float(LAYOUT-1)-0.5)
*
* ... Inside / Outside function defined by:
      fun=min(s1(xp,yp,zp),s2(xp,yp,zp),s3(xp,yp,zp),s4(xp,yp,zp),
     ,        s5(xp,yp,zp),s6(xp,yp,zp),s7(xp,yp,zp),s8(xp,yp,zp))
* ... The function is: negative for points outside the cube,
* .................... positive for points inside the cube,
* .................... zero for points at the cube's faces.
*
* ... Uniform field:
      write(*,*) fun
20    continue
15    continue
10    continue
      end