This subroutine is locating the facets of elements that are on the patch described by a 4-9 node quadrilateral surface. The routine works by computing a rectangular bounding box to the surface. The box is set based on a planar surface (z(1,*) and z(2,*)) and then computing the maxiumum and minimum value for this plane (z(3,*)) values. The programmer chose to use a hierarchical form for the interpolations to do the Newton iteration to find the values -- this is a choice and does not have anything to do with theory. The remainder of the routine then searches for facets that are within the bounding box. This routine could fail if the surface search has values for min/max z(3,*) too large (praticularly in polar/cylindrical form) but for small variations is quite reliable. Currently, also, it only works for 4-node tets and 8-node bricks.
