* function or manually set neighbor pointer from and to those
   * elements.
   */
  void unset_has_neighbor_ptrs()
  { _preparation.has_neighbor_ptrs = false; }

  /**
   * Tells this we have done some operation (e.g. adding elements with

   * User code which adds new elements to this mesh must call this
   * function.
   */
  void unset_has_cached_elem_data()
  { _preparation.has_cached_elem_data = false; }

  /**
   * Tells this we have done some operation (e.g. refining elements

   * User code which moves nodes ... should probably call this method,
   * in case ghosting functors depending on position exist?
   */
  void unset_has_reinit_ghosting_functors()
  { _preparation.has_reinit_ghosting_functors = false; }

  /**
   * Tells this we have done some operation which may have invalidated

   * User code which removes elements, or which adds or removes
   * boundary entries, should call this method.
   */
  void unset_has_boundary_id_sets()
  { _preparation.has_boundary_id_sets = false; }

  /**
   * \returns \p true if all elements and nodes of the mesh

                                bool store_parent_side_ids,
                                const std::vector<subdomain_id_type> & new_subdomain_ids)
{
  libmesh_assert(new_subdomain_ids.empty() ||
                 new_subdomain_ids.size() == 1 ||
                 new_subdomain_ids.size() == requested_boundary_ids.size());
  LOG_SCOPE("add_elements()", "BoundaryInfo");

      for (auto s : elem->side_index_range())
        {
          bool add_this_side = false;
          boundary_id_type triggering_bcid = invalid_id;

          // Find all the boundary side ids for this Elem side.
          std::vector<boundary_id_type> bcids;

              if (requested_boundary_ids.count(bcid))
                {
                  add_this_side = true;
                  triggering_bcid = bcid;
                  break;
                }
            }

            add_this_side = true;

          if (add_this_side)
            sides_to_add.emplace_back(elem->id(), s, triggering_bcid);
        }
    }

  // contiguous vector lets us do that index lookup in O(1) (after an
  // O(log N) Utility::binary_find()) rather than the O(N) std::distance()
  // between std::set iterators.
  std::vector<boundary_id_type> requested_boundary_ids_vec;
  if (new_subdomain_ids.size() > 1)
    requested_boundary_ids_vec.assign(requested_boundary_ids.begin(),
                                      requested_boundary_ids.end());

  for (const auto & [elem_id, s, triggering_bcid] : sides_to_add)
    {
      Elem * elem = _mesh->elem_ptr(elem_id);

        new_elem->set_extra_integer(parent_side_index_tag, s);

      // Assign subdomain id if requested.
      if (new_subdomain_ids.size() == 1)
        new_elem->subdomain_id() = new_subdomain_ids[0];
      else if (new_subdomain_ids.size() > 1)
        {
          auto it = Utility::binary_find(requested_boundary_ids_vec.begin(),
                                         requested_boundary_ids_vec.end(),
                                         triggering_bcid);
          libmesh_assert(it != requested_boundary_ids_vec.end());
          new_elem->subdomain_id() =
            new_subdomain_ids[std::distance(requested_boundary_ids_vec.begin(), it)];
        }

#ifdef LIBMESH_ENABLE_AMR

#endif

  // global containers are not synced
  boundary_mesh.unset_has_boundary_id_sets();

  // New elements were added; mark the affected preparation flags invalid.
  boundary_mesh.unset_has_neighbor_ptrs();
  boundary_mesh.unset_has_cached_elem_data();
  boundary_mesh.unset_has_reinit_ghosting_functors();
}

    {
      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