// this captures the case that the point is not (almost) in the direction of the element.:
  // first, project the problem onto the unit sphere:
  Point p_o(p - my_origin);
  pt0_o /= pt0_o.norm();
  pt1_o /= pt1_o.norm();
  pt2_o /= pt2_o.norm();
  p_o /= p_o.norm();

  // now, check if it is in the projected face; by comparing the distance of
  // any point in the element to \p p with the largest distance between this point
  // to any other point in the element.
  if ((p_o - pt0_o).norm_sq() > std::max((pt0_o - pt1_o).norm_sq(), (pt0_o - pt2_o).norm_sq()) ||
      (p_o - pt1_o).norm_sq() > std::max((pt1_o - pt2_o).norm_sq(), (pt1_o - pt0_o).norm_sq()) ||
      (p_o - pt2_o).norm_sq() > std::max((pt2_o - pt0_o).norm_sq(), (pt2_o - pt1_o).norm_sq()) )
    {
      // the physical point is definitely not contained in the element
      return false;
    }

      _communicator.set_union(_es_id_to_name);
    }

  this->synchronize_global_id_set();
}



void BoundaryInfo::synchronize_global_id_set()
{
  // Handle global data
  _global_boundary_ids = _boundary_ids;
  libmesh_assert(_mesh);
  if (!_mesh->is_serial())
    _communicator.set_union(_global_boundary_ids);
}

void BoundaryInfo::remove_edge (const Elem * elem,
                                const unsigned short int edge,
                                const boundary_id_type id)
{
  libmesh_assert(elem);

  // Only touch BCs for edges that exist.
  libmesh_assert_less (edge, elem->n_edges());

  // Only level 0 elements are stored in BoundaryInfo.
  libmesh_assert_equal_to (elem->level(), 0);

  // Erase (elem, edge, id) entries from map.
  erase_if(_boundary_edge_id, elem,
           [edge, id](decltype(_boundary_edge_id)::mapped_type & pr)
           {return pr.first == edge && pr.second == id;});
}


void BoundaryInfo::remove_shellface (const Elem * elem,

{
  // Pass global==false to the sub-methods here, so we can avoid
  // unnecessary communications
  this->remove_side_id(id, false);
  this->remove_edge_id(id, false);
  this->remove_shellface_id(id, false);
  this->remove_node_id(id, false);

  if (global)
    _global_boundary_ids.erase(id);
}



void BoundaryInfo::remove_side_id (boundary_id_type id, const bool global)
{
  // Erase id from ids containers
  _side_boundary_ids.erase(id);

  if (!_edge_boundary_ids.count(id) &&
      !_shellface_boundary_ids.count(id) &&
      !_node_boundary_ids.count(id))
    _boundary_ids.erase(id);

  _ss_id_to_name.erase(id);
  if (global)
    {
      bool someone_has_it = _boundary_ids.count(id);
      this->comm().max(someone_has_it);
      if (!someone_has_it)
        _global_boundary_ids.erase(id);
    }

  // Erase (*, *, id) entries from map.
  erase_if(_boundary_side_id,
           [id](decltype(_boundary_side_id)::mapped_type & pr)
           {return pr.second == id;});
}



void BoundaryInfo::remove_edge_id (boundary_id_type id, const bool global)
{
  // Erase id from ids containers
  _edge_boundary_ids.erase(id);

  if (!_side_boundary_ids.count(id) &&
      !_shellface_boundary_ids.count(id) &&
      !_node_boundary_ids.count(id))
    _boundary_ids.erase(id);

  _es_id_to_name.erase(id);
  if (global)
    {
      bool someone_has_it = _boundary_ids.count(id);
      this->comm().max(someone_has_it);
      if (!someone_has_it)
        _global_boundary_ids.erase(id);
    }

  // Erase (*, *, id) entries from map.
  erase_if(_boundary_edge_id,
           [id](decltype(_boundary_edge_id)::mapped_type & pr)
           {return pr.second == id;});
}



void BoundaryInfo::remove_shellface_id (boundary_id_type id, const bool global)
{
  // Erase id from ids containers
  _shellface_boundary_ids.erase(id);

  if (!_side_boundary_ids.count(id) &&
      !_edge_boundary_ids.count(id) &&
      !_node_boundary_ids.count(id))
    _boundary_ids.erase(id);

  if (global)
    {
      bool someone_has_it = _boundary_ids.count(id);
      this->comm().max(someone_has_it);
      if (!someone_has_it)
        _global_boundary_ids.erase(id);
    }

  // Erase (*, *, id) entries from map.
  erase_if(_boundary_shellface_id,
           [id](decltype(_boundary_shellface_id)::mapped_type & pr)
           {return pr.second == id;});
}



void BoundaryInfo::remove_node_id (boundary_id_type id, const bool global)
{
  // Erase id from ids containers
  _node_boundary_ids.erase(id);

  if (!_side_boundary_ids.count(id) &&
      !_edge_boundary_ids.count(id) &&
      !_shellface_boundary_ids.count(id))
    _boundary_ids.erase(id);

  if (global)
    {
      bool someone_has_it = _boundary_ids.count(id);
      this->comm().max(someone_has_it);
      if (!someone_has_it)
        _global_boundary_ids.erase(id);
    }

  // Erase (*, id) entries from map.
  erase_if(_boundary_node_id,
           [id](decltype(_boundary_node_id)::mapped_type & val)
           {return val == id;});
}

  renumber_name(_ns_id_to_name, old_id, new_id);
  renumber_name(_es_id_to_name, old_id, new_id);

  this->libmesh_assert_valid_multimaps();
}



void BoundaryInfo::renumber_side_id (boundary_id_type old_id,
                                     boundary_id_type new_id)
{
  // If the IDs are the same, this is a no-op.
  if (old_id == new_id)
    return;

  bool found_side = false;
  for (auto & p : _boundary_side_id)
    if (p.second.second == old_id)
      {
        // If we already have this id on this side, we don't want to
        // create a duplicate in our multimap
        this->remove_side(p.first, p.second.first, new_id);
        p.second.second = new_id;
        found_side = true;
      }
  _side_boundary_ids.erase(old_id);

  if (found_side)
    {
      _side_boundary_ids.insert(new_id);

      if (!_shellface_boundary_ids.count(old_id) &&
          !_edge_boundary_ids.count(old_id) &&
          !_node_boundary_ids.count(old_id))
        {
          _boundary_ids.erase(old_id);
        }
      _boundary_ids.insert(new_id);
    }

  renumber_name(_ss_id_to_name, old_id, new_id);

  this->libmesh_assert_valid_multimaps();
}



void BoundaryInfo::renumber_edge_id (boundary_id_type old_id,
                                     boundary_id_type new_id)
{
  // If the IDs are the same, this is a no-op.
  if (old_id == new_id)
    return;

  bool found_edge = false;
  for (auto & p : _boundary_edge_id)
    if (p.second.second == old_id)
      {
        // If we already have this id on this edge, we don't want to
        // create a duplicate in our multimap
        this->remove_edge(p.first, p.second.first, new_id);
        p.second.second = new_id;
        found_edge = true;
      }
  _edge_boundary_ids.erase(old_id);

  if (found_edge)
    {
      _edge_boundary_ids.insert(new_id);

      if (!_shellface_boundary_ids.count(old_id) &&
          !_side_boundary_ids.count(old_id) &&
          !_node_boundary_ids.count(old_id))
        {
          _boundary_ids.erase(old_id);
        }
      _boundary_ids.insert(new_id);
    }

  renumber_name(_es_id_to_name, old_id, new_id);

  this->libmesh_assert_valid_multimaps();
}



void BoundaryInfo::renumber_shellface_id (boundary_id_type old_id,
                                          boundary_id_type new_id)
{
  // If the IDs are the same, this is a no-op.
  if (old_id == new_id)
    return;

  bool found_shellface = false;
  for (auto & p : _boundary_shellface_id)
    if (p.second.second == old_id)
      {
        // If we already have this id on this shellface, we don't want
        // to create a duplicate in our multimap
        this->remove_shellface(p.first, p.second.first, new_id);
        p.second.second = new_id;
        found_shellface = true;
      }
  _shellface_boundary_ids.erase(old_id);

  if (found_shellface)
    {
      _shellface_boundary_ids.insert(new_id);

      if (!_edge_boundary_ids.count(old_id) &&
          !_side_boundary_ids.count(old_id) &&
          !_node_boundary_ids.count(old_id))
        {
          _boundary_ids.erase(old_id);
        }
      _boundary_ids.insert(new_id);
    }

  this->libmesh_assert_valid_multimaps();
}



void BoundaryInfo::renumber_node_id (boundary_id_type old_id,
                                     boundary_id_type new_id)
{
  // If the IDs are the same, this is a no-op.
  if (old_id == new_id)
    return;

  bool found_node = false;
  for (auto & p : _boundary_node_id)
    if (p.second == old_id)
      {
        // If we already have this id on this node, we don't want to
        // create a duplicate in our multimap
        this->remove_node(p.first, new_id);
        p.second = new_id;
        found_node = true;
      }
  _node_boundary_ids.erase(old_id);

  if (found_node)
    {
      _node_boundary_ids.insert(new_id);

      if (!_shellface_boundary_ids.count(old_id) &&
          !_side_boundary_ids.count(old_id) &&
          !_edge_boundary_ids.count(old_id))
        {
          _boundary_ids.erase(old_id);
        }
      _boundary_ids.insert(new_id);
    }

  renumber_name(_ns_id_to_name, old_id, new_id);

  this->libmesh_assert_valid_multimaps();
}

void
BoundaryInfo::build_node_list_from_side_list(const std::set<boundary_id_type> & sideset_list)
{
  // If we're on a distributed mesh, even the owner of a node is not
  // guaranteed to be able to properly assign its new boundary id(s)!

    {
      // Don't add remote sides
      if (elem->is_remote())
        continue;

      auto [sidenum, bcid] = id_pair;

      if (!sideset_list.empty() && !sideset_list.count(bcid))
        continue;

      // Need to loop over the sides of any possible children
      std::vector<const Elem *> family;
#ifdef LIBMESH_ENABLE_AMR
      elem->active_family_tree_by_side (family, sidenum);
#else
      family.push_back(elem);
#endif

      for (const auto & cur_elem : family)
        {
          side = &side_builder(*cur_elem, sidenum);

          // Add each node node on the side with the side's boundary id
          for (auto i : side->node_index_range())

       node_id_action_functor, datum_type_ex);
}

void BoundaryInfo::build_side_list_from_node_list(const std::set<boundary_id_type> & nodeset_list)
{
  // Check for early return
  if (_boundary_node_id.empty())
    {
      libMesh::out << "No boundary node IDs have been added: cannot build side list!" << std::endl;
      return;
    }

  // For avoiding extraneous element side construction
  ElemSideBuilder side_builder;
  // Pull objects out of the loop to reduce heap operations
  const Elem * side_elem;

  for (const auto & elem : _mesh->active_element_ptr_range())
    for (auto side : elem->side_index_range())
      {
        side_elem = &side_builder(*elem, side);

        // map from nodeset_id to count for that ID
        std::map<boundary_id_type, unsigned> nodesets_node_count;

        // For each nodeset that this node is a member of, increment the associated
        // nodeset ID count
        for (const auto & node : side_elem->node_ref_range())
          for (const auto & pr : as_range(_boundary_node_id.equal_range(&node)))
            if (nodeset_list.empty() || nodeset_list.count(pr.second))
              nodesets_node_count[pr.second]++;

        // Now check to see what nodeset_counts have the correct
        // number of nodes in them.  For any that do, add this side to
        // the sideset, making sure the sideset inherits the
        // nodeset's name, if there is one.
        for (const auto & pr : nodesets_node_count)
          if (pr.second == side_elem->n_nodes())
            {
              add_side(elem, side, pr.first);

              // Let the sideset inherit any non-empty name from the nodeset
              std::string & nset_name = nodeset_name(pr.first);

              if (nset_name != "")
                sideset_name(pr.first) = nset_name;
            }
      } // end for side
}