std::unique_ptr<MeshBase>
BoundaryElementConversionGenerator::generate()
{
  UnstructuredMesh & mesh = dynamic_cast<UnstructuredMesh &>(*_input);

  libMesh::MeshSerializer serial_mesh(mesh);

  // collect the subdomain ids of the original mesh
  std::set<subdomain_id_type> original_subdomain_ids;

  }

  // MeshSerializer's destructor calls delete_remote_elements()
  if (!mesh.is_replicated())
    mesh.prepare_for_use();
  else
    mesh.set_isnt_prepared();

  return std::move(_input);
}

            *other_mesh, _positions[i](0), _positions[i](1), _positions[i](2));
      }

      MooseMeshUtils::copyIntoMesh(*this,
                                   *mesh,
                                   *other_mesh,
                                   _avoid_merging_subdomains,
                                   _avoid_merging_boundaries,
                                   _communicator);
    }

          *translated_mesh, _positions[i](0), _positions[i](1), _positions[i](2));

      // Copy into final mesh
      MooseMeshUtils::copyIntoMesh(*this,
                                   *final_mesh,
                                   *translated_mesh,
                                   _avoid_merging_subdomains,
                                   _avoid_merging_boundaries,
                                   _communicator);

      // Reset nodal coordinates

                        "Whether to convert the 3D hole meshes with non-TRI3 surface sides into "
                        "a compatible form.");

  MooseEnum conversion_method("ALL SURFACE", "ALL");
  params.addParam<MooseEnum>(
      "conversion_method",
      conversion_method,
      "The method to convert 3D hole meshes into compatible meshes. Options are "
      "ALL: convert all elements into TET4; SURFACE: convert only the surface elements "
      "that are stitched to the generated Delaunay mesh.");

  params.addParam<bool>(
      "combined_stitching",
      false,
      "Whether to stitch all holes in one combined stitching step. This is efficient if a great "
      "number of holes are to be stitched. But it is hard to debug if problems arise.");

    _hole_ptrs(getMeshes("holes")),
    _stitch_holes(getParam<std::vector<bool>>("stitch_holes")),
    _convert_holes_for_stitching(getParam<bool>("convert_holes_for_stitching")),
    _conversion_method(parameters.get<MooseEnum>("conversion_method")),
    _combined_stitching(parameters.get<bool>("combined_stitching")),
    _algorithm(parameters.get<MooseEnum>("algorithm")),
    _verbose_stitching(parameters.get<bool>("verbose_stitching"))
{

      paramError("hole_boundaries", "Need one hole_boundaries entry per hole, if specified.");
  }

  if (isParamSetByUser("conversion_method") && !_convert_holes_for_stitching)
    paramError(
        "conversion_method",
        "This parameter is only applicable when convert_holes_for_stitching is set to true.");
}

    // We do not need to worry about element order here as libMesh checks it
    std::set<ElemType> hole_elem_types;
    std::set<unsigned short> hole_elem_dims;
    std::vector<std::pair<dof_id_type, unsigned int>> hole_elem_external_sides;
    for (auto elem : hole_mesh.element_ptr_range())
    {
      hole_elem_dims.emplace(elem->dim());

          if (!elem->neighbor_ptr(s))
          {
            hole_elem_types.emplace(elem->side_ptr(s)->type());
            hole_elem_external_sides.emplace_back(elem->id(), s);
          }
        }
      // For a non-3D element, we just need to record the element type

          paramError("holes",
                     "3D hole meshes with non-TRI3 surface elements cannot be stitched without "
                     "converting them to TET4. Consider setting convert_holes_for_stitching=true.");
        else if (_stitch_holes.size() && _stitch_holes[hole_i] && _conversion_method == "SURFACE")
        {
          // Create a transition layer with triangle sides on the external boundary of the hole
          BoundaryID temp_ext_bid = MooseMeshUtils::getNextFreeBoundaryID(hole_mesh);
          for (const auto & hees : hole_elem_external_sides)
            hole_mesh.get_boundary_info().add_side(hees.first, hees.second, temp_ext_bid);
          MooseMeshElementConversionUtils::transitionLayerGenerator(
              hole_mesh, std::vector<BoundaryName>({std::to_string(temp_ext_bid)}), 1, false);
          hole_mesh.get_boundary_info().remove_id(temp_ext_bid);
          // MeshSerializer's destructor calls delete_remote_elements()
          // For replicated meshes, we need to refresh the neighbor information
          if (!hole_mesh.is_replicated())
            hole_mesh.prepare_for_use();
          else
            hole_mesh.find_neighbors();
        }
        else
          MeshTools::Modification::all_tri(**_hole_ptrs[hole_i]);

      const auto & increment_subdomain_map = hole_mesh.get_subdomain_name_map();
      main_subdomain_map.insert(increment_subdomain_map.begin(), increment_subdomain_map.end());

      if (_combined_stitching)
        MooseMeshUtils::copyIntoMesh(*this, *mesh, hole_mesh, false, false, _communicator);
      else
      {
        std::size_t n_nodes_stitched = mesh->stitch_meshes(hole_mesh,
                                                           inner_bcid,
                                                           new_hole_bcid,
                                                           TOLERANCE,
                                                           /*clear_stitched_bcids*/ true,
                                                           _verbose_stitching,
                                                           use_binary_search);

        if (!n_nodes_stitched)
          mooseError("Failed to stitch hole mesh ", hole_i, " to new tetrahedralization.");
      }
    }
  }
  if (doing_stitching && _combined_stitching)
  {
    std::size_t n_nodes_stitched = mesh->stitch_surfaces(inner_bcid,
                                                         new_hole_bcid,
                                                         TOLERANCE,
                                                         /*clear_stitched_bcids*/ true,
                                                         _verbose_stitching,
                                                         use_binary_search);
    if (!n_nodes_stitched)
      mooseError("Failed to stitch combined hole meshes to new tetrahedralization.");
  }

  // Add user-specified sideset names

namespace MooseMeshElementConversionUtils
{
void
hexElemSplitter(MeshBase & mesh,
                const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
                const dof_id_type elem_id,
                std::vector<dof_id_type> & converted_elems_ids)

}

void
prismElemSplitter(MeshBase & mesh,
                  const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
                  const dof_id_type elem_id,
                  std::vector<dof_id_type> & converted_elems_ids)

}

void
pyramidElemSplitter(MeshBase & mesh,
                    const std::vector<libMesh::BoundaryInfo::BCTuple> & bdry_side_list,
                    const dof_id_type elem_id,
                    std::vector<dof_id_type> & converted_elems_ids)

}

void
convert3DMeshToAllTet4(MeshBase & mesh,
                       const std::vector<std::pair<dof_id_type, bool>> & elems_to_process,
                       std::vector<dof_id_type> & converted_elems_ids_to_track,
                       const subdomain_id_type block_id_to_remove,

}

void
convert3DMeshToAllTet4(MeshBase & mesh)
{
  mooseAssert(mesh.is_serial(),
              "This method only supports serial meshes. If you are calling this method with a "

}

void
convertElem(MeshBase & mesh,
            const dof_id_type & elem_id,
            const std::vector<unsigned int> & side_indices,
            const std::vector<std::vector<boundary_id_type>> & elem_side_info,

}

void
convertHex8Elem(MeshBase & mesh,
                const dof_id_type & elem_id,
                const std::vector<unsigned int> & side_indices,
                const std::vector<std::vector<boundary_id_type>> & elem_side_info,

}

void
createUnitPyramid5FromHex8(MeshBase & mesh,
                           const dof_id_type & elem_id,
                           const unsigned int & side_index,
                           const Node * new_node,

}

void
createUnitTet4FromHex8(MeshBase & mesh,
                       const dof_id_type & elem_id,
                       const unsigned int & side_index,
                       const Node * new_node,

}

void
convertPrism6Elem(MeshBase & mesh,
                  const dof_id_type & elem_id,
                  const std::vector<unsigned int> & side_indices,
                  const std::vector<std::vector<boundary_id_type>> & elem_side_info,

}

void
createUnitTet4FromPrism6(MeshBase & mesh,
                         const dof_id_type & elem_id,
                         const unsigned int & side_index,
                         const Node * new_node,

}

void
createUnitPyramid5FromPrism6(MeshBase & mesh,
                             const dof_id_type & elem_id,
                             const unsigned int & side_index,
                             const Node * new_node,

}

void
convertPyramid5Elem(MeshBase & mesh,
                    const dof_id_type & elem_id,
                    const std::vector<std::vector<boundary_id_type>> & elem_side_info,
                    const SubdomainID & subdomain_id_shift_base)

}

void
retainEEID(MeshBase & mesh, const dof_id_type & elem_id, Elem * new_elem_ptr)
{
  const unsigned int n_eeid = mesh.n_elem_integers();
  for (const auto & i : make_range(n_eeid))

}

void
transitionLayerGenerator(MeshBase & mesh,
                         const std::vector<BoundaryName> & boundary_names,
                         const unsigned int conversion_element_layer_number,
                         const bool external_boundaries_checking)

}

void
copyIntoMesh(MeshGenerator & mg,
             UnstructuredMesh & destination,
             const UnstructuredMesh & source,
             const bool avoid_merging_subdomains,
             const bool avoid_merging_boundaries,
             const Parallel::Communicator & communicator)
{
  dof_id_type node_delta = destination.max_node_id();
  dof_id_type elem_delta = destination.max_elem_id();

  unique_id_type unique_delta =
#ifdef LIBMESH_ENABLE_UNIQUE_ID
      destination.parallel_max_unique_id();
#else
      0;
#endif

  // Prevent overlaps by offsetting the subdomains in
  std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
  unsigned int block_offset = 0;
  if (avoid_merging_subdomains)
  {
    // Note: if performance becomes an issue, this is overkill for just getting the max node id
    std::set<subdomain_id_type> source_ids;
    std::set<subdomain_id_type> dest_ids;
    source.subdomain_ids(source_ids, true);
    destination.subdomain_ids(dest_ids, true);
    mooseAssert(source_ids.size(), "Should have a subdomain");
    mooseAssert(dest_ids.size(), "Should have a subdomain");
    unsigned int max_dest_bid = *dest_ids.rbegin();
    unsigned int min_source_bid = *source_ids.begin();
    communicator.max(max_dest_bid);
    communicator.min(min_source_bid);
    block_offset = 1 + max_dest_bid - min_source_bid;
    for (const auto bid : source_ids)
      id_remapping[bid] = block_offset + bid;
  }

  // Copy mesh data over from the other mesh
  destination.copy_nodes_and_elements(source,
                                      // Skipping this should cause the neighbors
                                      // to simply be copied from the other mesh
                                      // (which makes sense and is way faster)

                                      avoid_merging_subdomains ? &id_remapping : nullptr);

  // Get an offset to prevent overlaps / wild merging between boundaries
  BoundaryInfo & boundary = destination.get_boundary_info();
  const BoundaryInfo & other_boundary = source.get_boundary_info();

  unsigned int bid_offset = 0;
  if (avoid_merging_boundaries)
  {
    const auto boundary_ids = boundary.get_boundary_ids();
    const auto other_boundary_ids = other_boundary.get_boundary_ids();
    unsigned int max_dest_bid = boundary_ids.size() ? *boundary_ids.rbegin() : 0;
    unsigned int min_source_bid = other_boundary_ids.size() ? *other_boundary_ids.begin() : 0;
    communicator.max(max_dest_bid);
    communicator.min(min_source_bid);
    bid_offset = 1 + max_dest_bid - min_source_bid;
  }

  // Note: the code below originally came from ReplicatedMesh::stitch_mesh_helper()
  // in libMesh replicated_mesh.C around line 1203

  // Copy BoundaryInfo from other_mesh too.  We do this via the
  // list APIs rather than element-by-element for speed.
  for (const auto & t : other_boundary.build_node_list())
    boundary.add_node(std::get<0>(t) + node_delta, bid_offset + std::get<1>(t));

  for (const auto & t : other_boundary.build_side_list())
    boundary.add_side(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));

  for (const auto & t : other_boundary.build_edge_list())
    boundary.add_edge(std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));

  for (const auto & t : other_boundary.build_shellface_list())
    boundary.add_shellface(
        std::get<0>(t) + elem_delta, std::get<1>(t), bid_offset + std::get<2>(t));

  // Check for the case with two block ids sharing the same name
  if (avoid_merging_subdomains)
  {
    mooseAssert(mg.parameters().isParamDefined("avoid_merging_subdomains"),
                "Missing parameter in the mesh generator calling this function: "
                "avoid_merging_subdomains. Considering setting avoid_merging_subdomains to true.");
    for (const auto & [block_id, block_name] : destination.get_subdomain_name_map())
      for (const auto & [source_id, source_name] : source.get_subdomain_name_map())
        if (block_name == source_name)
          mg.paramWarning(
              "avoid_merging_subdomains",
              "Not merging subdomains is creating two subdomains with the same name '" +
                  block_name + "' but different ids: " + std::to_string(source_id) + " & " +
                  std::to_string(block_id + block_offset) +
                  ".\n We recommend using a RenameBlockGenerator to prevent this as you "
                  "will get errors reading the Exodus output later.");
  }

  for (const auto & [block_id, block_name] : source.get_subdomain_name_map())
    destination.set_subdomain_name_map().insert(
        std::make_pair<SubdomainID, SubdomainName>(block_id + block_offset, block_name));

  // Check for the case with two boundary ids sharing the same name
  if (avoid_merging_boundaries)
  {
    mooseAssert(mg.parameters().isParamDefined("avoid_merging_boundaries"),
                "Missing parameter in the mesh generator calling this function: "
                "avoid_merging_boundaries. Considering setting avoid_merging_boundaries to true.");
    for (const auto & [b_id, b_name] : other_boundary.get_sideset_name_map())
      for (const auto & [source_id, source_name] : boundary.get_sideset_name_map())
        if (b_name == source_name)
          mg.paramWarning(
              "avoid_merging_boundaries",
              "Not merging boundaries is creating two sidesets with the same name '" + b_name +
                  "' but different ids: " + std::to_string(source_id) + " & " +
                  std::to_string(b_id + bid_offset) +
                  ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
                  "will get errors reading the Exodus output later.");
    for (const auto & [b_id, b_name] : other_boundary.get_nodeset_name_map())
      for (const auto & [source_id, source_name] : boundary.get_nodeset_name_map())
        if (b_name == source_name)
          mg.paramWarning(
              "avoid_merging_boundaries",
              "Not merging boundaries is creating two nodesets with the same name '" + b_name +
                  "' but different ids: " + std::to_string(source_id) + " & " +
                  std::to_string(b_id + bid_offset) +
                  ".\n We recommend using a RenameBoundaryGenerator to prevent this as you "
                  "will get errors reading the Exodus output later.");
  }

  for (const auto & [nodeset_id, nodeset_name] : other_boundary.get_nodeset_name_map())
    boundary.set_nodeset_name_map().insert(
        std::make_pair<BoundaryID, BoundaryName>(nodeset_id + bid_offset, nodeset_name));

  for (const auto & [sideset_id, sideset_name] : other_boundary.get_sideset_name_map())
    boundary.set_sideset_name_map().insert(
        std::make_pair<BoundaryID, BoundaryName>(sideset_id + bid_offset, sideset_name));

  for (const auto & [edgeset_id, edgeset_name] : other_boundary.get_edgeset_name_map())
    boundary.set_edgeset_name_map().insert(
        std::make_pair<BoundaryID, BoundaryName>(edgeset_id + bid_offset, edgeset_name));
}
}