/**
   * Const getter for the _point attribute
   */
  const Point &point() const { return _point; }

  /**
   * Const getter for the _tol attribute

 * Represents an invalid constraint (i.e., when the two constraints don't
 * intersect)
 */
class InvalidConstraint
{

public:
  InvalidConstraint()
    : _err_msg("We should never get here! The InvalidConstraint object should be "
               "detected and replaced with a valid ConstraintVariant prior to calling "
               "any class methods.")
  {
  }

  /**
   * Dummy intersect method that should never be called.

 * @param b Constraint to combine with a.
 * @return Combination (intersection) of constraint a and b.
 */
inline ConstraintVariant intersect_constraints(const ConstraintVariant &a,
                                               const ConstraintVariant &b) {
  // std::visit applies the visitor v (a Callable that can be called with any
  // combination of types from Variants) to the active value inside a

  // ConstraintVariant do)
  return std::visit(
      [](const auto &lhs, const auto &rhs) -> ConstraintVariant {
        return lhs.intersect(rhs);
      },
      a, b);
}

                      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);

auto get_positive_vector = [](const Point & vec) -> Point {
  // unit() rounds to within a few ulps of length 1, so TOLERANCE only
  // rejects vectors that were never normalized
  libmesh_assert(absolute_fuzzy_equals(vec.norm(), 1., TOLERANCE));

  // Choose sign such that direction vector points in the positive x/y/z
  // direction This helps to eliminate duplicate lines/planes

  const auto dot_prod = vec * canonical;
  libmesh_assert(!absolute_fuzzy_equals(dot_prod, 0.));

  return (dot_prod > 0) ? vec : Point(-vec);
};

PointConstraint::PointConstraint(const Point & point, const Real & tol) : _point(point), _tol(tol)

LineConstraint::LineConstraint(const Point & point, const Point & direction, const Real & tol)
  : _point(point), _direction(get_positive_vector(direction)), _tol(tol)
{
  libmesh_assert_msg(_direction.norm() >= _tol,
                     "Can't define a line with zero magnitude direction vector.");
}

bool LineConstraint::is_parallel(const PlaneConstraint & p) const
{
  return std::abs(_direction * p.normal()) < _tol;
}

ConstraintVariant LineConstraint::intersect(const ConstraintVariant & other) const

PlaneConstraint::PlaneConstraint(const Point & point, const Point & normal, const Real & tol)
  : _point(point), _normal(get_positive_vector(normal)), _tol(tol)
{
  libmesh_assert_msg(_normal.norm() >= _tol,
                     "Can't define a plane with zero magnitude direction vector.");
}

  // |n1 x n2| is the sine of the angle between the unit normals, and its
  // square is the determinant that intersect() divides by, so any pair of
  // planes that fails this test can be intersected.
  return _normal.cross(p.normal()).norm() <= _tol;
}

bool PlaneConstraint::is_parallel(const LineConstraint & l) const { return l.is_parallel(*this); }

            // nonzero because the planes are not parallel. Computing it from
            // dir avoids the cancellation in 1 - cos^2 for nearly parallel
            // planes.
            const Real denom = -dir.norm_sq();

            const Real s = -(n1_dot_n2 * n2_dot_w - n2_dot_n2 * n1_dot_w) / denom;
            const Point p0 = _point + s * _normal;

            return LineConstraint{p0, dir.unit()};
          }

        else if constexpr (std::is_same_v<T, LineConstraint>)

            if (this->is_parallel(o))
              // Line is parallel and does not intersect the plane
              return InvalidConstraint();

            // Solve for t in the parametric equation:
            //   p(t) = point + t·d

              if (candidate_normal.norm() <= TOLERANCE)
                continue;

              const PlaneConstraint candidate_plane{node, candidate_normal.unit()};
              valid_planes.emplace(candidate_plane);
            }
        }

      // Fall back to fixed node constraint
      if (std::holds_alternative<InvalidConstraint>(current))
        {
          current = PointConstraint(node);
          break;
        }
    }