Fix bounding ball calculation for coaxial lumped ports

[sebastiangrimberg]

Previously, for coaxial lumped ports with particularly coarse meshes, the following assertion would fail in geodata.cpp's BoundingBallFromPointCloud:

MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <=
                rel_tol * ball.radius,
            "Furthest point perpendicular must be on the exterior of the ball: "
                << PerpendicularDistance({delta}, perp) << " vs. " << ball.radius
                << "!");

This PR improves the bounding ball algorithm to avoid this error (and is a simpler algorithm). We now just compute the ball origin as the average of all points and radial direction as the point furthest from the centroid. The rest of the algorithm is unchanged.

Here is an example second-order mesh which was causing problems prior to this PR:

[sebastiangrimberg]

Converted to draft while revising the bounding box/bounding ball algorithms.

[sebastiangrimberg]

Commit c22bf04 corrects the previous commit after noticing that the oriented bounding box calculation from mesh::GetBoundingBox only provides a correct result when the convex hull of the points forms a rectangle or rectangular prism. For spheres/circles, the box is inscribed in the circle/sphere so the use in CoaxialElementData needs adjusting.

This could be an opportunity to implement a true oriented bounding box calculation, if it might be useful elsewhere. One option is from Baraquet and Har-Peled, where the core ingredient would be a 2D minimal bounding rectangle implementation. See, for example, the implementation here which is Boost Licensed.

Another idea is to first compute the convex hull of the points, followed by an SVD-based approach for computing the approximate oriented bounding box.

Alternatively, for now, the current implementation works fine for boundaries which are rectangles and circles as we expect in current usage.

[sebastiangrimberg]

Following up on this comment, the latest PR implements Welzl's algorithm to compute the minimum bounding sphere which is used for coaxial ports now. The comment for when the simplified OBB calculation given by mesh::ComputeBoundingBox has been clarified.

[hughcars]

LGTM. Some minor notes on possibly swapping, and to document that there isn't really an issue with the case where p_1 p_2 are p_3 and p_4.

[hughcars]

Note for posterity: By dropping the projection onto the chord in the search, p_3 and p_4 might be p_1 and p_2. This would cause issues with the algorithm for four points on a sphere. However, the Welzl function would only find two entries, then search the remaining space finding no more.

In practice this shouldn't matter, but could be fixed either by using projected distance, or excluding the p_1 and p_2 from the search space. When p_1 and p_2 are found, rather than add to the end, swap them to the end, and search over the reduced space of .end() - 2.

      std::iter_swap(p_1, vertices.end());
      std::iter_swap(p_2, vertices.end() - 1);
      auto p_3 =
          std::max_element(vertices.begin(), vertices.end()-2,
                           [&p_12](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
                           { return (x - p_12).norm() < (y - p_12).norm(); });
      auto p_4 = std::max_element(vertices.begin(), vertices.end()-2,
                                  [p_3](const Eigen::Vector3d &x, const Eigen::Vector3d &y)
                                  { return (x - *p_3).norm() < (y - *p_3).norm(); });

should suffice to ensure that p_3 and p_4 are different and found earlier.

[sebastiangrimberg]

I'll just revert to the projected distance, since it's already implemented. Not difficult to do and should be an improvement even if it doesn't matter much if at all for performance.

[sebastiangrimberg]

Done in 32ebd0c.

[hughcars]

Alternative: size exactly then iter_swap, p_i to the end and std::iota(vertices.begin(), vertices.end(), indices.begin(), 0);

[sebastiangrimberg]

Done in 32ebd0c.