program quad
*             ====
      parameter (LAYOUT=31)
* ... Fortran statement functions,
* ... defining the shape functions of a
* ... Finite Element 4-node quadrilateral:
      s1(x,y)=(0.5-x)*(0.5-y)
      s2(x,y)=(0.5+x)*(0.5-y)
      s3(x,y)=(0.5-x)*(0.5+y)
      s4(x,y)=(0.5+x)*(0.5+y)
*
* ... 31 x 31 x 31 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)
*
* ... Inside / Outside function defined by:
      fun=min(s1(xp,yp),s2(xp,yp),s3(xp,yp),s4(xp,yp))
* ... The function is: negative for points outside the quad,
* .................... positive for points inside the quad,
* .................... zero for points at the quad's faces.
*
* ... Uniform field:
      write(*,*) fun
15    continue
10    continue
      end