`Intersection2DPolygon` can return clockwise overlap vertices for nearly coincident triangles

[ebchin]

Running the contact patch example in Smith uncovered fragility in determining CCW orientation in Intersection2DPolygon() for the general intersection path. When intersecting polygons are nearly overlapping, but not overlapping enough to trigger the all A in B or all B in A paths, then PolyReorderConvex() can fail. Nearly overlapping polygons have some vertices that are nearly identical, which aren't merged until a subsequent call to CheckPolySegs(). When calling PolyReorderConvex() with nearly overlapping vertices, the vertices potentially do not lie on the boundary of a convex polygon, breaking the assumptions of PolyReorderConvex().

Issue reproducer:

#include "tribol/geom/GeomUtilities.hpp"

#include <iostream>
#include <iomanip>

static double signedArea(const double* x, const double* y, int n)
{
  double a = 0.0;
  for (int i = 0; i < n; ++i) {
    const int j = (i + 1) % n;
    a += x[i] * y[j] - y[i] * x[j];
  }
  return 0.5 * a;
}

int main()
{
  // Nonmortar projected triangle from failing contact_patch case.
  const double xA[3] = {
      -1.7720926303008128e-01,
       1.7720926303008130e-01,
       7.6327832942979648e-17
  };
  const double yA[3] = {
      -5.9069754343360363e-02,
      -5.9069754343360349e-02,
       1.1813950868672080e-01
  };

  // Mortar projected triangle after ElemReverse() in ComputeMortarForceEnzyme.
  const double xB[3] = {
       1.7720926303008075e-01,
      -8.3266726846886605e-17,
      -1.7720926303008128e-01
  };
  const double yB[3] = {
      -5.9069754343361050e-02,
       1.1813950868672063e-01,
      -5.9069754343360363e-02
  };

  double polyX[8] = {};
  double polyY[8] = {};
  int numPolyVert = 0;
  double area = 0.0;

  auto err = tribol::Intersection2DPolygon(
      xA, yA, 3,
      xB, yB, 3,
      1.0e-8, 1.0e-8,
      polyX, polyY,
      numPolyVert,
      area,
      true);

  std::cout << std::scientific << std::setprecision(16);
  std::cout << "err = " << err << "\n";
  std::cout << "numPolyVert = " << numPolyVert << "\n";
  std::cout << "area = " << area << "\n";
  std::cout << "signed_area = " << signedArea(polyX, polyY, numPolyVert) << "\n";
  std::cout << "poly = ";
  for (int i = 0; i < numPolyVert; ++i) {
    std::cout << (i == 0 ? "" : "; ") << polyX[i] << ", " << polyY[i];
  }
  std::cout << "\n";

  return signedArea(polyX, polyY, numPolyVert) > 0.0 ? 0 : 1;
}

And here is the output (we expect the output to be a CCW-oriented polygon):

err = 0
numPolyVert = 3
area = 3.1403122903664482e-02
signed_area = -3.1403122903664482e-02
orientation = CW
poly = -1.7720926303008128e-01, -5.9069754343360363e-02; -8.3266726846886605e-17, 1.1813950868672063e-01; 1.7720926303008130e-01, -5.9069754343360349e-02