}

inline const std::pair<unsigned int, unsigned int> &
DofMap::get_variable_array(const unsigned int vi) const
{
  auto it = std::upper_bound(
      _array_variables.begin(),
      _array_variables.end(),
      vi,
      [](unsigned int value, const std::pair<unsigned int, unsigned int> & b) { return value < b.first; });

  libmesh_assert_msg(it != _array_variables.begin(),
                     "Passed in " << std::to_string(vi) << " is not in any of our array variables");
  --it;
  libmesh_assert_msg(vi < it->second,
                     "Passed in " << std::to_string(vi) << " is not in any of our array variables");
  return *it;
}

template <typename DofIndicesFunctor>
void DofMap::array_dof_indices(const DofIndicesFunctor & functor,
                               std::vector<dof_id_type> & di,
                               const unsigned int vn) const
{
  const auto [begin, end] = this->get_variable_array(vn);
  functor(di, begin);

  const unsigned int count = end - begin;
  // We make count, which could be >> ntest, the inner index in hopes of vectorization
  if (count > 1)
    {
      const dof_id_type component_size = di.size();
      di.resize(count * component_size);

      const auto pack_container = [&di,
                                   component_size](const unsigned int j,
                                                   const std::vector<dof_id_type> & j_dof_indices,
                                                   const unsigned int stride) {
        if (&j_dof_indices != &di)
          libmesh_assert(j_dof_indices.size() == component_size);
        for (const auto i : make_range(component_size))
          di[j * component_size + i] = j_dof_indices[i] + stride * j;
      };
      pack_container(0, di, 0);

      const auto & fe_type = _variable_groups[libmesh_map_find(_var_to_vg, vn)].type();
      if (const bool lagrange = fe_type.family == LAGRANGE;
          lagrange || (FEInterface::get_continuity(fe_type) == DISCONTINUOUS))
        {
          const auto stride = lagrange ? 1 : component_size;
          for (const auto j : make_range((unsigned int)1, count))
            pack_container(j, di, stride);
        }
      else
        {
          static thread_local std::vector<dof_id_type> work_dof_indices;
          unsigned int j = 1;
          for (const auto i : make_range(begin + 1, end))
            {
              functor(work_dof_indices, i);
              pack_container(j++, work_dof_indices, 0);
            }
        }
    }
}

inline
unsigned int DofMap::n_vars() const
{
  return cast_int<unsigned int>(_variables.size());
}

inline
const std::string & DofMap::variable_name (const unsigned int i) const
{
  libmesh_assert_less (i, _variables.size());

  return _variables[i].name();
}

inline
bool DofMap::identify_variable_groups () const
{
  return _identify_variable_groups;
}

inline
void DofMap::identify_variable_groups (const bool ivg)
{
  _identify_variable_groups = ivg;
}

inline
unsigned int DofMap::n_components(const MeshBase & mesh) const
{
  if (_variables.empty())
    return 0;

  const Variable & last = _variables.back();
  return last.first_scalar_number() + last.n_components(mesh);
}

inline
unsigned int
DofMap::variable_scalar_number (unsigned int var_num,
                                unsigned int component) const
{
  return _variables[var_num].first_scalar_number() + component;
}

inline
const FEType & DofMap::variable_type (std::string_view var) const
{
  return _variables[this->variable_number(var)].type();
}

inline bool DofMap::has_variable(std::string_view var) const
{
  return _variable_numbers.count(var);
}

inline unsigned int DofMap::variable_number(std::string_view var) const
{
  auto var_num = libmesh_map_find(_variable_numbers, var);
  libmesh_assert_equal_to(_variables[var_num].name(), var);
  return var_num;
}

} // namespace libMesh

  _variable_groups.clear();
  _var_to_vg.clear();
  _variable_group_numbers.clear();
  _array_variables.clear();
  _first_scalar_df.clear();
  this->clear_send_list();
  this->clear_sparsity();

    bool all_discontinuous_dofs = true;

    // We may call this method even without ever having initialized our data
      for (auto var : index_range(this->_variables))
        if (FEInterface::get_continuity(this->variable_type(var)) !=  DISCONTINUOUS)
          all_discontinuous_dofs = false;

    if (all_discontinuous_dofs)
      implicit_neighbor_dofs = true;

         std::vector<dof_id_type> & dof_indices,
         const dof_id_type dof) { dof_indices.push_back(dof); },
      p_level);
}

void DofMap::array_dof_indices(const Elem * const elem,
                               std::vector<dof_id_type> & di,
                               const unsigned int vn,
                               int p_level) const
{
  auto dof_indices_functor = [elem, p_level, this](std::vector<dof_id_type> & functor_di,
                                                   const unsigned int functor_vn) {
    this->dof_indices(elem, functor_di, functor_vn, p_level);
  };
  this->array_dof_indices(dof_indices_functor, di, vn);
}

void DofMap::array_dof_indices(const Node * const node,
                               std::vector<dof_id_type> & di,
                               const unsigned int vn) const
{
  auto dof_indices_functor = [node, this](std::vector<dof_id_type> & functor_di,
                                          const unsigned int functor_vn) {
    this->dof_indices(node, functor_di, functor_vn);
  };
  this->array_dof_indices(dof_indices_functor, di, vn);
}

void DofMap::dof_indices (const Node * const node,
                          std::vector<dof_id_type> & di) const

    _sc->reinit();
}

unsigned int DofMap::add_variable(System & sys,
                                  std::string_view var,
                                  const FEType & type,
                                  const std::set<subdomain_id_type> * const active_subdomains)
{
  parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync

  libmesh_assert(this->comm().verify(std::string(var)));
  libmesh_assert(this->comm().verify(type));
  libmesh_assert(this->comm().verify((active_subdomains == nullptr)));

  if (active_subdomains)
    libmesh_assert(this->comm().verify(active_subdomains->size()));

  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  for (auto v : make_range(this->n_vars()))
    if (this->variable_name(v) == var)
      {
        if (this->variable_type(v) == type)
          {
            // Check whether the existing variable's active subdomains also matches
            // the incoming variable's active subdomains. If they don't match, then
            // either it is an error by the user or the user is trying to change the
            // subdomain restriction after the variable has already been added, which
            // is not supported.
            const Variable & existing_var = this->variable(v);

            // Check whether active_subdomains is not provided/empty and the existing_var is
            // implicitly_active()
            bool check1 = (!active_subdomains || active_subdomains->empty()) &&
                          existing_var.implicitly_active();

            // Check if the provided active_subdomains is equal to the existing_var's
            // active_subdomains
            bool check2 =
                (active_subdomains && (*active_subdomains == existing_var.active_subdomains()));

            // If either of these checks passed, then we already have this variable
            if (check1 || check2)
              return _variables[v].number();
          }

        libmesh_error_msg("ERROR: incompatible variable "
                          << var << " has already been added for this system!");
      }

  libmesh_assert(!sys.is_initialized());

  if (this->n_variable_groups())
    {
      // Optimize for VariableGroups here - if the user is adding multiple
      // variables of the same FEType and subdomain restriction, catch

      //
      // start by setting this flag to whatever the user has requested
      // and then consider the conditions which should negate it.
      bool should_be_in_vg = this->identify_variable_groups();

      VariableGroup & vg = _variable_groups.back();

      // get a pointer to their subdomain restriction, if any.
      const std::set<subdomain_id_type> * const their_active_subdomains(
          vg.implicitly_active() ? nullptr : &vg.active_subdomains());

      // Different types?
      if (vg.type() != type)
        should_be_in_vg = false;

      // they are restricted, we aren't?
      if (their_active_subdomains &&
          (!active_subdomains || (active_subdomains && active_subdomains->empty())))
        should_be_in_vg = false;

      // they aren't restricted, we are?
      if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
        should_be_in_vg = false;

      if (their_active_subdomains && active_subdomains)
        // restricted to different sets?
        if (*their_active_subdomains != *active_subdomains)
          should_be_in_vg = false;

      // OK, after all that, append the variable to the vg if none of the conditions
      // were violated
      if (should_be_in_vg)
        {
          const unsigned int vn = this->n_vars();

          std::string varstr(var);

          _variable_numbers[varstr] = vn;
          vg.append(std::move(varstr));
          _variables.push_back(vg(vg.n_variables() - 1));
          const unsigned int vgn = _variable_groups.size() - 1;
          _variable_group_numbers.push_back(vgn);
          _var_to_vg.emplace(vn, vgn);

          return vn;
        }
    }

  // otherwise, fall back to adding a single variable group
  return this->add_variables(
      sys, std::vector<std::string>(1, std::string(var)), type, active_subdomains);
}

unsigned int DofMap::add_variables(System & sys,
                                   const std::vector<std::string> & vars,
                                   const FEType & type,
                                   const std::set<subdomain_id_type> * const active_subdomains)
{
  parallel_object_only(); // Not strictly needed, but the only safe way to keep in sync

  libmesh_assert(!sys.is_initialized());

  libmesh_assert(this->comm().verify(vars.size()));
  libmesh_assert(this->comm().verify(type));
  libmesh_assert(this->comm().verify((active_subdomains == nullptr)));

  if (active_subdomains)
    libmesh_assert(this->comm().verify(active_subdomains->size()));

  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  for (auto ovar : vars)
    {
      libmesh_assert(this->comm().verify(ovar));

      for (auto v : make_range(this->n_vars()))
        if (this->variable_name(v) == ovar)
          {
            if (this->variable_type(v) == type)
              return _variables[v].number();

            libmesh_error_msg("ERROR: incompatible variable "
                              << ovar << " has already been added for this system!");
          }
    }

  if (this->n_variable_groups())
    {
      // Optimize for VariableGroups here - if the user is adding multiple
      // variables of the same FEType and subdomain restriction, catch

      //
      // start by setting this flag to whatever the user has requested
      // and then consider the conditions which should negate it.
      bool should_be_in_vg = this->identify_variable_groups();

      VariableGroup & vg = _variable_groups.back();

      // get a pointer to their subdomain restriction, if any.
      const std::set<subdomain_id_type> * const their_active_subdomains(
          vg.implicitly_active() ? nullptr : &vg.active_subdomains());

      // Different types?
      if (vg.type() != type)
        should_be_in_vg = false;

      // they are restricted, we aren't?
      if (their_active_subdomains &&
          (!active_subdomains || (active_subdomains && active_subdomains->empty())))
        should_be_in_vg = false;

      // they aren't restricted, we are?
      if (!their_active_subdomains && (active_subdomains && !active_subdomains->empty()))
        should_be_in_vg = false;

      if (their_active_subdomains && active_subdomains)
        // restricted to different sets?
        if (*their_active_subdomains != *active_subdomains)
          should_be_in_vg = false;

      // If after all that none of the conditions were violated,
      // append the variables to the vg and we're done
      if (should_be_in_vg)
        {
          unsigned int vn = this->n_vars();
          const unsigned int vgn = _variable_groups.size() - 1;

          for (auto ovar : vars)
            {
              vn = this->n_vars();

              vg.append(ovar);

              _variables.push_back(vg(vg.n_variables() - 1));
              _variable_numbers[ovar] = vn;
              _variable_group_numbers.push_back(vgn);
              _var_to_vg.emplace(vn, vgn);
            }
          return vn;
        }
    }

  const unsigned int curr_n_vars = this->n_vars();

  const unsigned int next_first_component = this->n_components(sys.get_mesh());

  // We weren't able to add to an existing variable group, so
  // add a new variable group to the list
  _variable_groups.push_back(
      (active_subdomains == nullptr)
          ? VariableGroup(&sys, vars, curr_n_vars, next_first_component, type)
          : VariableGroup(&sys, vars, curr_n_vars, next_first_component, type, *active_subdomains));

  const VariableGroup & vg(_variable_groups.back());
  const unsigned int vgn = _variable_groups.size() - 1;

  // Add each component of the group individually
  for (auto v : make_range(vars.size()))
    {
      const unsigned int vn = curr_n_vars + v;
      _variables.push_back(vg(v));
      _variable_numbers[vars[v]] = vn;
      _variable_group_numbers.push_back(vgn);
      _var_to_vg.emplace(vn, vgn);
    }

  libmesh_assert_equal_to((curr_n_vars + vars.size()), this->n_vars());

  // BSK - Defer this now to System::init_data() so we can detect
  // VariableGroups 12/28/2012

  // _dof_map->add_variable_group (vg);

  // Return the number of the new variable
  return cast_int<unsigned int>(curr_n_vars + vars.size() - 1);
}

unsigned int DofMap::add_variable_array (System & sys,
                                         const std::vector<std::string> & vars,
                                         const FEType & type,
                                         const std::set<subdomain_id_type> * const active_subdomains)
{
  const unsigned int count = cast_int<unsigned int>(vars.size());
  const unsigned int last_var = this->add_variables(sys, vars, type, active_subdomains);
  const unsigned int first_var = last_var + 1 - count;
  _array_variables.push_back({first_var, first_var + count});
  return last_var;
}

void DofMap::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
{
  all_variable_numbers.resize(n_vars());

  unsigned int count = 0;
  for (auto vn : _variable_numbers)
    all_variable_numbers[count++] = vn.second;
}

template LIBMESH_EXPORT bool DofMap::is_evaluable<Elem>(const Elem &, unsigned int) const;
template LIBMESH_EXPORT bool DofMap::is_evaluable<Node>(const Node &, unsigned int) const;

                              if (neighbor->p_level() < my_p_level &&
                                  neighbor->p_refinement_flag() != Elem::REFINE)
                                {
                                  neighbor->set_p_refinement_flag(Elem::REFINE);
                                  level_one_satisfied = false;
                                  compatible_with_coarsening = false;
                                }
                              if (neighbor->p_level() == my_p_level &&

    {
      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

                                   const FEType & type,
                                   const std::set<subdomain_id_type> * const active_subdomains)
{
  return this->get_dof_map().add_variable(*this, var, type, active_subdomains);
}

                                    const FEType & type,
                                    const std::set<subdomain_id_type> * const active_subdomains)
{
  return this->get_dof_map().add_variables(*this, vars, type, active_subdomains);
}

                             active_subdomains);
}

unsigned int System::add_variable_array (const std::vector<std::string> & vars,
                                         const FEType & type,
                                         const std::set<subdomain_id_type> * const active_subdomains)
{
  return this->get_dof_map().add_variable_array(*this, vars, type, active_subdomains);
}

bool System::has_variable (std::string_view var) const
{
  return this->get_dof_map().has_variable(var);
}

unsigned int System::variable_number (std::string_view var) const
{
  return this->get_dof_map().variable_number(var);
}

void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
{
  this->get_dof_map().get_all_variable_numbers(all_variable_numbers);
}

  return this->get_dof_map().has_static_condensation();
}

unsigned int System::n_vars() const
{
  return this->get_dof_map().n_vars();
}

const std::string & System::variable_name (const unsigned int i) const
{
  return this->get_dof_map().variable_name(i);
}

bool System::identify_variable_groups () const
{
  return this->get_dof_map().identify_variable_groups();
}

void System::identify_variable_groups (const bool ivg)
{
  this->get_dof_map().identify_variable_groups(ivg);
}

unsigned int System::n_components() const
{
  return this->get_dof_map().n_components(this->get_mesh());
}

unsigned int System::n_variable_groups() const
{
  return this->get_dof_map().n_variable_groups();
}

const Variable & System::variable (const unsigned int i) const
{
  return this->get_dof_map().variable(i);
}

const VariableGroup & System::variable_group (const unsigned int vg) const
{
  return this->get_dof_map().variable_group(vg);
}

unsigned int
System::variable_scalar_number (unsigned int var_num,
                                unsigned int component) const
{
  return this->get_dof_map().variable_scalar_number(var_num, component);
}

const FEType & System::variable_type (const unsigned int i) const
{
  return this->get_dof_map().variable_type(i);
}

const FEType & System::variable_type (std::string_view var) const
{
  return this->get_dof_map().variable_type(var);
}

} // namespace libMesh