};


inline void ExodusII_IO_Helper::set_add_sides(bool add_sides)
{
  _add_sides = add_sides;
}


inline bool ExodusII_IO_Helper::get_add_sides()

  if (old_n)
    {
      *old_n = p;
      old_n->processor_id() = proc_id;

      return old_n;
    }

  Node * n = Node::build(p, id).release();

  {
    // Without a reduction across ranks, the subdomain map is only
    // complete if the whole mesh is available here.
    libmesh_assert(reduce || mesh.is_serial_on_zero());
    if (reduce)
      libmesh_parallel_only(mesh.comm());

    // When reducing, each rank only contributes its local elements, so
    // that no element is counted twice on a replicated mesh.
    auto for_each_elem = [&mesh, reduce](auto && f)
      {
        if (reduce)
          for (const auto & elem : mesh.active_local_element_ptr_range())
            f(*elem);
        else
          for (const auto & elem : mesh.active_element_ptr_range())
            f(*elem);
      };

    // Exodus requires every element block to contain a single element
    // type.  Group elements first by subdomain id and then by element

    // (for example a mix of C0POLYHEDRON and HEX8 cells) is split across
    // multiple blocks rather than written as a single, invalid block.
    std::map<subdomain_id_type,
             std::map<ElemType, std::vector<unsigned int>>> elems_by_subdomain_type;

    // We skip writing infinite elements to the Exodus file, so don't
    // put them in the map. That way the number of blocks should be
    // correct.
    if (reduce)
      {
        std::vector<dof_id_type> elem_ids;
        std::vector<subdomain_id_type> elem_sbd_ids;
        std::vector<int> elem_types;
        for_each_elem([&](const Elem & elem)
          {
            if (elem.infinite())
              return;
            elem_ids.push_back(elem.id());
            elem_sbd_ids.push_back(elem.subdomain_id());
            elem_types.push_back(elem.type());
          });

        // Every rank contributes the same number of entries to each
        // vector, so the gathered vectors stay aligned.
        mesh.comm().allgather(elem_ids);
        mesh.comm().allgather(elem_sbd_ids);
        mesh.comm().allgather(elem_types);

        for (auto i : index_range(elem_ids))
          elems_by_subdomain_type[elem_sbd_ids[i]][static_cast<ElemType>(elem_types[i])]
            .push_back(cast_int<unsigned int>(elem_ids[i]));

        // Gathered ids arrive grouped by rank; sort them to match the
        // id ordering of a serial traversal.
        for (auto & [sbd_id, type_map] : elems_by_subdomain_type)
          {
            libmesh_ignore(sbd_id);
            for (auto & [elem_t, ids] : type_map)
              {
                libmesh_ignore(elem_t);
                std::sort(ids.begin(), ids.end());
              }
          }
      }
    else
      for_each_elem([&](const Elem & elem)
        {
          if (elem.infinite())
            return;
          elems_by_subdomain_type[elem.subdomain_id()][elem.type()].push_back(elem.id());
        });

    // Assign a block id to each (subdomain, element type) group.  The

    // subdomain get synthesized block ids allocated above all existing
    // subdomain ids.
    std::map<subdomain_id_type, std::vector<unsigned int>> subdomain_map;
    subdomain_id_type next_synth_block_id = elems_by_subdomain_type.empty() ?
      0 : elems_by_subdomain_type.rbegin()->first + 1;

    for (auto & [sbd_id, type_map] : elems_by_subdomain_type)
      {
        bool first_type = true;
        for (auto & [elem_t, elem_ids] : type_map)
          {
            libmesh_ignore(elem_t);
            const subdomain_id_type block_id =
              first_type ? sbd_id : next_synth_block_id++;

            first_type = false;
            subdomain_map[block_id] = std::move(elem_ids);
          }
      }

    // Real element blocks occupy the ids below subdomain_id_end; any
    // blocks synthesized for visualization sides are numbered after them.
    if (!subdomain_map.empty())
      subdomain_id_end = subdomain_map.rbegin()->first + 1;

    // If we've been asked to add side elements, those go in their own
    // blocks.  We don't have any ids to list for elements that don't

    // the number of elements we'll add in each new block.
    if (add_sides)
      {
        std::map<subdomain_id_type, unsigned int> side_block_counts;
        for_each_elem([&](const Elem & elem)
          {
            if (elem.infinite())
              return;

            for (auto s : elem.side_index_range())
              {
                if (EquationSystems::redundant_added_side(elem,s))
                  continue;

                const subdomain_id_type side_block_id =
                  cast_int<subdomain_id_type>(subdomain_id_end + elem.side_type(s));

                // Guard the invariant above (also catches subdomain_id_type
                // overflow wrapping a side block id back down onto an
                // element block id).
                libmesh_assert_greater_equal(side_block_id, subdomain_id_end);

                ++side_block_counts[side_block_id];
              }
          });

        if (reduce)
          mesh.comm().sum(side_block_counts);

        for (const auto & [side_block_id, count] : side_block_counts)
          subdomain_map[side_block_id] = {count};
      }

    return subdomain_map;

    return;

  // We need to know about all processors' subdomains
  subdomain_id_type subdomain_id_end = 0;
  auto subdomain_map = build_subdomain_map(mesh,
                                           _add_sides,
                                           /*reduce=*/!_run_only_on_proc0,
                                           subdomain_id_end);

  if (!use_discontinuous)
    {

{
  LOG_SCOPE("write_elements()", "ExodusII_IO_Helper");

  if ((_run_only_on_proc0) && (this->processor_id() != 0))
    return;

  // Map from block ID to a vector of element IDs in that block.  Element
  // IDs are now of type dof_id_type, subdomain IDs are of type subdomain_id_type.
  subdomain_id_type subdomain_id_end = 0;
  auto subdomain_map = build_subdomain_map(mesh,
                                           _add_sides,
                                           /*reduce=*/!_run_only_on_proc0,
                                           subdomain_id_end);

  // Allocate the polyhedron face block ID after both mesh subdomains and
  // any blocks synthesized for visualization sides.
  const int c0polyhedron_face_block_id = subdomain_map.empty() ? 0 :
    cast_int<int>(subdomain_map.rbegin()->first) + 1;

  // element map vector
  num_elem_blk = cast_int<int>(subdomain_map.size());

          // detect such synthesized blocks from the mismatch between the
          // block id and the actual subdomain of the elements it holds.
          const subdomain_id_type elem_subdomain_id =
            mesh.elem_ref(element_id_vec[0]).subdomain_id();
          const bool is_synthesized_block = (subdomain_id != elem_subdomain_id);

          std::string block_name = mesh.subdomain_name(elem_subdomain_id);

          if (is_synthesized_block)
            {
              // Suffix the (original) subdomain's name with the element
              // type, so the synthesized block is identifiable and its
              // name stays distinct from the block that kept the
              // subdomain id.  This block becomes its own subdomain when
              // the file is read back in.
              const std::string type_suffix = Utility::enum_to_string<ElemType>(elem_t);
              block_name = block_name.empty()
                ? (std::to_string(elem_subdomain_id) + "_" + type_suffix)
                : (block_name + "_" + type_suffix);

              // Informational: keep it to one rank so parallel (Nemesis)
              // writes will not repeat it once per processor.
              // NOTE this might miss logs if only occurs on other ranks
              if (this->processor_id() == 0)
                libMesh::out << "ExodusII_IO: subdomain " << elem_subdomain_id
                             << " contains more than one element type; writing its "
                             << type_suffix << " elements to a new block \"" << block_name
                             << "\" (block id " << subdomain_id
                             << ")"
                             << std::endl;
            }

          if (block_name.empty() && elem_t == C0POLYGON)

                  continue;

                if (elem->side_type(s) != elem_t)
                  continue;

                const std::vector<unsigned int> side_nodes =
                  elem->nodes_on_side(s);

                                  const bool synchronous)
{
  if (synchronous)
    parallel_object_only();
  else
    this->unset_has_synched_subdomain_name_map();
  _block_id_to_name[id] =  name;

    {
      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

  // if the node already exists, then assign new (x,y,z) values
  if (n)
    *n = p;
  // otherwise build a new node, put it in the right spot, and return
  // a valid pointer.
  else

#ifndef NDEBUG
      if (elem->runtime_topology() && elem->dim() == 3)
        for (auto s : elem->side_index_range())
          libmesh_assert(elem->nodes_on_side(s) ==
                         polyhedron_side_nodes[s]);
#endif