Surface integration tools appear to be coded to first order elements

Look at: https://github.com/sandialabs/optimism/blob/4f94d76d88c7a16b17f4fff4f1f51fb5ed922151/optimism/Surface.py

I think it assumes linear elements in several functions.

Fails silently for high-order elements. Contact iterations do not converge, but no reason is given.