* 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.

  Point v = this->point(1) - this->point(0);
  if (!v.relative_fuzzy_equals
      ((this->point(4) - this->point(0))*2, affine_tol))
    return false;
  v = this->point(2) - this->point(1);
  if (!v.relative_fuzzy_equals
      ((this->point(5) - this->point(1))*2, affine_tol))
    return false;
  v = this->point(2) - this->point(0);
  if (!v.relative_fuzzy_equals
      ((this->point(6) - this->point(0))*2, affine_tol))

  v = this->point(3) - this->point(0);
  if (!v.relative_fuzzy_equals
      ((this->point(7) - this->point(0))*2, affine_tol))
    return false;
  v = this->point(3) - this->point(1);
  if (!v.relative_fuzzy_equals
      ((this->point(8) - this->point(1))*2, affine_tol))

  Point v10 = this->point(1) - this->point(0);
  if (!v10.relative_fuzzy_equals
      ((this->point(4) - this->point(0))*2))
    return false;
  Point v21 = this->point(2) - this->point(1);
  if (!v21.relative_fuzzy_equals
      ((this->point(5) - this->point(1))*2))

  find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
}

void find_nodal_or_face_neighbors(
    const MeshBase & mesh,
    const Node & node,
    const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,

{
  // Find all the nodal neighbors... that is the nodes directly connected
  // to this node through one edge.
  find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);

  // If no neighbors are found, use all nodes on the containing side as
  // neighbors.
  if (!neighbors.size())
    {
      // Grab the element containing node
      const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
      // Find the element side containing node
      for (const auto &side : elem->side_index_range())
        {
          const auto &nodes_on_side = elem->nodes_on_side(side);
          const auto it =
              std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
                return elem->node_id(local_node_id) == node.id();
              });

          if (it != nodes_on_side.end())
            {
              for (const auto &local_node_id : nodes_on_side)
                // No need to add node itself as a neighbor
                if (const auto *node_ptr = elem->node_ptr(local_node_id);
                    *node_ptr != node)
                  neighbors.push_back(node_ptr);
              break;
            }
        }
    }
  libmesh_assert(neighbors.size());
}

            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

  // to this node through one edge, or if none exists, use other nodes on the
  // containing face
  std::vector<const Node *> neighbors;
  MeshTools::find_nodal_or_face_neighbors(
      mesh, node, nodes_to_elem_map, neighbors);

  // Each constituent set corresponds to neighbors sharing a face on the

  // to this node through one edge, or if none exists, use other nodes on the
  // containing face
  std::vector<const Node *> neighbors;
  MeshTools::find_nodal_or_face_neighbors(
      mesh, node, nodes_to_elem_map, neighbors);

  // Each constituent set corresponds to neighbors sharing a face on the

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

  const auto ref_vol = target_elem->reference_elem()->volume();

  // Update the nodes of the target element, depending on type
  const Real sqrt_2 = std::sqrt(Real(2));
  const Real sqrt_3 = std::sqrt(Real(3));
  std::vector<std::unique_ptr<Node>> owned_nodes;

    } // if Pyramid

  // Elems deriving from Tet
  else if (type_str.compare(0, 3, "TET") == 0)
    {

      // The ideal target element is a a regular tet with equilateral

      // of the non-optimal reference element (i.e., a right tet).

      // Side length that preserves the volume of the reference element
      const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol);
      // tet height with the property that all faces are equilateral triangles
      const auto target_height = sqrt_2 / sqrt_3 * side_length;

      const auto & s = side_length;
      const auto & h = target_height;

      // For regular tet
      //                                         x        y                z     node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s,  h),  3));

      if (type == TET10 || type == TET14)
        {
          const auto & on = owned_nodes;
          // Define the edge midpoint nodes of the tet

          // Base node to base node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6));
          // Base node to apex node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9));

          if (type == TET14)
            {
              // Define the face midpoint nodes of the tet
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10));
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11));
              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12));
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13));
            }
        }

      else if (type != TET4)
        libmesh_error_msg("Unsupported tet element: " << type_str);

    } // if Tet