// which is Knupp's algebraic skew (1 is ideal).
      // See: C. J. Stimpson et al., "The Verdict Geometric Quality
      // Library," Sandia report SAND2007-1751, 2007.
    case SKEW_ANGLE:
      {
        const Point
          x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),
          x4 = point(4), x5 = point(5), x6 = point(6), x7 = point(7);

        // Principal axes, one per logical (xi, eta, zeta) direction.
        const Point
          X1 = (x1 - x0) + (x2 - x3) + (x5 - x4) + (x6 - x7),
          X2 = (x3 - x0) + (x2 - x1) + (x7 - x4) + (x6 - x5),
          X3 = (x4 - x0) + (x5 - x1) + (x6 - x2) + (x7 - x3);

        const Real n1 = X1.norm(), n2 = X2.norm(), n3 = X3.norm();

        // Degenerate element: return 0 (the Verdict convention) if any
        // principal axis has zero length.
        if (n1 == 0. || n2 == 0. || n3 == 0.)
          return 0.;

        // Normalize, then take the largest |cos| among the three
        // pairs of principal axes.
        const Point X1h = X1 / n1, X2h = X2 / n2, X3h = X3 / n3;
        return std::max({std::abs(X1h * X2h),
                         std::abs(X1h * X3h),
                         std::abs(X2h * X3h)});
      }
#endif // LIBMESH_DIM >= 3

      bounds.second = 4.;
      break;

    case SKEW_ANGLE:
      bounds.first  = 0.;
      bounds.second = 0.5;
      break;

    case SKEW:
    case SHEAR:
    case SHAPE:
      bounds.first  = 0.3;

Real Elem::length(const unsigned int n1,
                  const unsigned int n2) const
{
  libmesh_assert_less ( n1, this->n_vertices() );
  libmesh_assert_less ( n2, this->n_vertices() );

  return (this->point(n1) - this->point(n2)).norm();
}

      // available to this per-element method, so we use the reference
      // element instead. Unlike the standard Verdict metric, J is not
      // squared here.
    case SIZE:
      {
        // 1D elements don't have interior corners, so this metric does
        // not really apply to them.
        const auto N = this->dim();
        if (N < 2)
          return 1.;

        // Average the nodal Jacobian determinant over the corner
        // nodes. This uses the same nodal Jacobian construction as the
        // JACOBIAN metric above.
        Real sum_node_area = 0.;
        unsigned int n_corners = 0;

        for (auto n : this->node_index_range())
          {
            // Get list of edge ids adjacent to this node.
            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);

            // Skip any nodes that don't have dim() adjacent edges (see
            // the JACOBIAN metric above for the Pyramid apex caveat).
            if (adjacent_edge_ids.size() != N)
              continue;

            // Construct oriented edges pointing away from node n.
            std::vector<Point> oriented_edges(N);
            for (auto i : make_range(N))
              {
                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
                if (node_0 != n)
                  std::swap(node_0, node_1);
                oriented_edges[i] = this->point(node_1) - this->point(node_0);
              }

            // Unscaled nodal area (2D) or volume (3D).
            Real node_area = (N == 2) ?
              cross_norm(oriented_edges[0], oriented_edges[1]) :
              std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));

            sum_node_area += node_area;
            ++n_corners;
          }

        // No usable corners, or a degenerate (zero-size) element: return
        // 0 (the lowest quality).
        if (n_corners == 0)
          return 0.;

        const Real J = sum_node_area / n_corners;
        if (J == 0.)
          return 0.;

        return std::min(J, Real(1) / J);
      }

      // Maximum condition number of the nodal Jacobian matrix over the

      // determinant) has an infinite condition number, reported as 0
      // following the convention used elsewhere (e.g. EDGE_LENGTH_RATIO)
      // that 0 stands in for infinity.
    case CONDITION:
      {
        // 1D elements don't have interior corners, so this metric does
        // not really apply to them.
        const auto N = this->dim();
        if (N < 2)
          return 1.;

        // kappa >= 1 for every matrix, so 1 is both the ideal value and
        // a safe floor for the running maximum.
        Real max_cond = 1.;

        for (auto n : this->node_index_range())
          {
            // Get list of edge ids adjacent to this node.
            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);

            // Skip any nodes that don't have dim() adjacent edges (see
            // the JACOBIAN metric above for the Pyramid apex caveat).
            if (adjacent_edge_ids.size() != N)
              continue;

            // Construct oriented edges pointing away from node n; these
            // are the columns of the nodal Jacobian A.
            std::vector<Point> e(N);
            for (auto i : make_range(N))
              {
                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
                if (node_0 != n)
                  std::swap(node_0, node_1);
                e[i] = this->point(node_1) - this->point(node_0);
              }

            // Squared Frobenius norm of A.
            Real frob_A_sq = 0.;
            for (auto i : make_range(N))
              frob_A_sq += e[i].norm_sq();

            // |det(A)| and the squared Frobenius norm of A^{-1}.
            Real abs_det, frob_Ainv_sq;
            if (N == 2)
              {
                abs_det = cross_norm(e[0], e[1]);

                // Degenerate corner: infinite condition number.
                if (abs_det == 0.)
                  return 0.;

                // For a 2x2 matrix, |A^{-1}|_F = |A|_F / |det|.
                frob_Ainv_sq = frob_A_sq / (abs_det * abs_det);
              }
            else
              {
                abs_det = std::abs(triple_product(e[0], e[1], e[2]));

                // Degenerate corner: infinite condition number.
                if (abs_det == 0.)
                  return 0.;

                // The rows of A^{-1} are (e1 x e2), (e2 x e0), (e0 x e1),
                // each divided by det(A).
                frob_Ainv_sq = (e[1].cross(e[2]).norm_sq() +
                                e[2].cross(e[0]).norm_sq() +
                                e[0].cross(e[1]).norm_sq()) / (abs_det * abs_det);
              }

            const Real kappa = std::sqrt(frob_A_sq * frob_Ainv_sq) / N;
            max_cond = std::max(max_cond, kappa);
          }

        return max_cond;
      }

      // Return 1 if we made it here

      its_name = "Skew";
      break;

    case SKEW_ANGLE:
      its_name = "Skew Angle";
      break;

    case SHEAR:
      its_name = "Shear";

           << '\n'
           << "Suggested ranges:\n"
           << "Hexes: (0.3 -> 1)\n"
           << "Quads: (0.3 -> 1)";
      break;

    case SKEW_ANGLE:
      desc << "Maximum |cos A|, where A\n"
           << "is the angle between edges\n"
           << "at element center.\n"

      // and return the smallest (worst) such ratio. The value lies in
      // (0, 1], with 1 indicating no taper, i.e. both pairs of opposite
      // edges are equal in length (as for any parallelogram).
    case TAPER:
      {
        const Real d01 = this->length(0,1);
        const Real d12 = this->length(1,2);
        const Real d23 = this->length(2,3);
        const Real d03 = this->length(0,3);

        // Longer length of each opposite-edge pair.
        const Real max0 = std::max(d01, d23);
        const Real max1 = std::max(d12, d03);

        // Degenerate element with a zero-length pair of opposite edges:
        // return 0 (the lowest quality).
        if (max0 == 0. || max1 == 0.)
          return 0.;

        return std::min(std::min(d01, d23) / max0,
                        std::min(d12, d03) / max1);
      }

    case SHAPE:

      // which is Knupp's algebraic skew (1 is ideal).
      // See: C. J. Stimpson et al., "The Verdict Geometric Quality
      // Library," Sandia report SAND2007-1751, 2007.
    case SKEW_ANGLE:
      {
        const Point x0 = this->point(0), x1 = this->point(1),
                    x2 = this->point(2), x3 = this->point(3);

        // Principal axes: midpoint-to-midpoint of opposite edges.
        const Point X1 = (x1 - x0) + (x2 - x3);
        const Point X2 = (x3 - x0) + (x2 - x1);

        const Real n1 = X1.norm(), n2 = X2.norm();

        // Degenerate element: return 0 (the Verdict convention) if
        // either principal axis has zero length.
        if (n1 == 0. || n2 == 0.)
          return 0.;

        return std::abs((X1 * X2) / (n1 * n2));
      }

      // This test returns 0 if a Quad:

      bounds.second = 1.;
      break;

    case SKEW_ANGLE:
      bounds.first  = 0.;
      bounds.second = 0.5;
      break;

    case DISTORTION:
      bounds.first  = 0.6;

                      sender_could_become_owner)
                    {
                      if (it != repartitioned_node_pids.end() &&
                          pid < it->second)
                        it->second = pid;
                      else
                        repartitioned_node_pids[n] = pid;
                    }
                  else
                    if (it == repartitioned_node_pids.end())
                      repartitioned_node_pids[n] =
                        DofObject::invalid_processor_id;

                  repartitioned_node_sets_to_push[pid].insert(n);

    {
      ray_target = inside - Point(1);
      intersection_distances =
        this->find_ray_intersections(inside, ray_target);
    }

  // I'd make this an assert, but I'm not 100% confident we can't