Draft: Nodal/Variationally Consistent Platen BCs

Extremely WIP! There is a lot of debugging code, logging, and hacked together solutions which will be removed later.

Currently, this MR implements variationally consistent platen BCs by directly looping over each nodal residual value and applying the nonlinear function

\bm{r} = \bm{r} - [\bm{r}\cdot (-\bm{n}_{platen}) - \gamma p(u)]_{\mathbb{R}^-} \bm{n}_{platen}

where \bm{n}_{platen} is the platen normal, p(u) is the penetration depth in the direction of the face normal \bm{n}, and \gamma is the Nitsche method parameter. The penetration p(u) is computed by first finding the penetration and conservative normals at each node using two surface QFunctions. The penetration normal is computed as

\bm{n}_p = \int_{S} -g(u) \gamma \phi(x) \hat{\bm{n}}(x)

and the conservative normal is computed as

\bar{\bm{n}} = \int_{S} \phi(x) \hat{\bm{n}}(x).

The surface area of each node is then A = ||\bar{\bm{n}}||_2, which is used to compute the (signed) penetration depth as \gamma p(u) = \bm{n}_p \cdot \bar{\bm{n}} / A. We compute the penetration depth in this way to allow for distinguishing between a gap, i.e. p(u)<0, and true penetration, i.e. p(u) > 0.

This formulation results in exact contact computations, seemingly regardless of the parameter \gamma. To verify, you can run any of the platen examples with the additional flags -multigrid none -snes_fd -fine_dm_mat_type seqaij.

I need help (@jedbrown and @jeremylt) figuring out the best way to change the static function RatelApplyPlatenBoundaryReactionForce in ratel-petsc-ops.c into a QFunction and nodal CeedOperator which uses the current residual and the updated nodal normal vectors. I know that we can only have one active input on the CeedOperator (which should probably be the residual), so is there another way that I can get the current normals into the operator, e.g. as passive inputs or otherwise?

P.S. I will go back and rebase out most of these commits/debug code later.