dof_map.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_DOF_MAP_H
#define LIBMESH_DOF_MAP_H
 
// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/libmesh.h" // libMesh::invalid_uint
#include "libmesh/variable.h"
#include "libmesh/threads.h"
#include "libmesh/threads_allocators.h"
#include "libmesh/elem_range.h"
#include "libmesh/ghosting_functor.h"
#include "libmesh/sparsity_pattern.h"
#include "libmesh/parallel_object.h"
#include "libmesh/point.h"
#include "libmesh/utility.h"
#include "libmesh/elem.h"
#include "libmesh/fe_interface.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/enum_elem_type.h"
#include "libmesh/mesh_subdivision_support.h"
#include "libmesh/dof_map_base.h"
 
// TIMPI includes
#include "timpi/parallel_implementation.h"
#include "timpi/parallel_sync.h"
 
// C++ Includes
#include <algorithm>
#include <cstddef>
#include <iterator>
#include <map>
#include <string>
#include <vector>
#include <memory>
 
namespace libMesh
{
 
// Forward Declarations
class CouplingMatrix;
class DefaultCoupling;
class DirichletBoundary;
class DirichletBoundaries;
class DofMap;
class DofObject;
class FEType;
class MeshBase;
class PeriodicBoundaryBase;
class PeriodicBoundaries;
class System;
class NonlinearImplicitSystem;
class StaticCondensationDofMap;
template <typename T> class DenseVectorBase;
template <typename T> class DenseVector;
template <typename T> class DenseMatrix;
template <typename T> class SparseMatrix;
template <typename T> class NumericVector;
enum Order : int;
 
 
 
// ------------------------------------------------------------
// Do we need constraints for anything?
 
#if defined(LIBMESH_ENABLE_AMR) ||              \
  defined(LIBMESH_ENABLE_PERIODIC) ||           \
  defined(LIBMESH_ENABLE_DIRICHLET)
#  define LIBMESH_ENABLE_CONSTRAINTS 1
#endif
 
// ------------------------------------------------------------
// AMR constraint matrix types
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
typedef std::map<dof_id_type, Real,
                 std::less<dof_id_type>,
                 Threads::scalable_allocator<std::pair<const dof_id_type, Real>>> DofConstraintRow;
 
class DofConstraints : public std::map<dof_id_type,
                                       DofConstraintRow,
                                       std::less<dof_id_type>,
                                       Threads::scalable_allocator<std::pair<const dof_id_type, DofConstraintRow>>>
{
};
 
class DofConstraintValueMap :
    public std::map<dof_id_type, Number,
                    std::less<dof_id_type>,
                    Threads::scalable_allocator<std::pair<const dof_id_type, Number>>>
{
};
 
class AdjointDofConstraintValues :
    public std::map<unsigned int, DofConstraintValueMap,
                    std::less<unsigned int>,
                    Threads::scalable_allocator
                    <std::pair<const unsigned int, DofConstraintValueMap>>>
{
};
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
typedef std::map<const Node *, Real,
                 std::less<const Node *>,
                 Threads::scalable_allocator<std::pair<const Node * const, Real>>> NodeConstraintRow;
 
class NodeConstraints : public std::map<const Node *,
                                        std::pair<NodeConstraintRow,Point>,
                                        std::less<const Node *>,
                                        Threads::scalable_allocator<std::pair<const Node * const, std::pair<NodeConstraintRow,Point>>>>
{
};
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
#endif // LIBMESH_ENABLE_CONSTRAINTS
 
 
 
class DofMap : public DofMapBase,
               public ReferenceCountedObject<DofMap>
{
public:
 
  explicit
  DofMap(const unsigned int sys_number,
         MeshBase & mesh);
 
  ~DofMap();
 
  class AugmentSparsityPattern : public SparsityPattern::AugmentSparsityPattern
  {};
 
  class AugmentSendList
  {
  public:
    virtual ~AugmentSendList () = default;
 
    virtual void augment_send_list (std::vector<dof_id_type> & send_list) = 0;
  };
 
  void attach_matrix (SparseMatrix<Number> & matrix);
 
  void update_sparsity_pattern(SparseMatrix<Number> & matrix) const;
 
  bool is_attached (SparseMatrix<Number> & matrix);
 
  std::size_t distribute_dofs (MeshBase &);
 
  void compute_sparsity (const MeshBase &);
 
  bool computed_sparsity_already () const;
 
  void set_constrained_sparsity_construction(bool use_constraints);
 
  void full_sparsity_pattern_needed();
 
  bool constrained_sparsity_construction();
 
  void clear_sparsity();
 
  void remove_default_ghosting();
 
  void add_default_ghosting();
 
  typedef std::vector<GhostingFunctor *>::const_iterator GhostingFunctorIterator;
 
  void add_coupling_functor(GhostingFunctor & coupling_functor,
                            bool to_mesh = true);
 
  void add_coupling_functor(std::shared_ptr<GhostingFunctor> coupling_functor,
                            bool to_mesh = true)
  { _shared_functors[coupling_functor.get()] = coupling_functor;
    this->add_coupling_functor(*coupling_functor, to_mesh); }
 
  void remove_coupling_functor(GhostingFunctor & coupling_functor);
 
  GhostingFunctorIterator coupling_functors_begin() const
  { return _coupling_functors.begin(); }
 
  GhostingFunctorIterator coupling_functors_end() const
  { return _coupling_functors.end(); }
 
  DefaultCoupling & default_coupling() { return *_default_coupling; }
 
  void add_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor,
                                      bool to_mesh = true);
 
  void add_algebraic_ghosting_functor(std::shared_ptr<GhostingFunctor> evaluable_functor,
                                      bool to_mesh = true)
  { _shared_functors[evaluable_functor.get()] = evaluable_functor;
    this->add_algebraic_ghosting_functor(*evaluable_functor, to_mesh); }
 
  void remove_algebraic_ghosting_functor(GhostingFunctor & evaluable_functor);
 
  GhostingFunctorIterator algebraic_ghosting_functors_begin() const
  { return _algebraic_ghosting_functors.begin(); }
 
  GhostingFunctorIterator algebraic_ghosting_functors_end() const
  { return _algebraic_ghosting_functors.end(); }
 
  DefaultCoupling & default_algebraic_ghosting() { return *_default_evaluating; }
 
  void attach_extra_sparsity_object (SparsityPattern::AugmentSparsityPattern & asp)
  {
    _augment_sparsity_pattern = &asp;
  }
 
  void attach_extra_sparsity_function(void (*func)(SparsityPattern::Graph & sparsity,
                                                   std::vector<dof_id_type> & n_nz,
                                                   std::vector<dof_id_type> & n_oz,
                                                   void *),
                                      void * context = nullptr)
  { _extra_sparsity_function = func; _extra_sparsity_context = context; }
 
  void attach_extra_send_list_object (DofMap::AugmentSendList & asl)
  {
    _augment_send_list = &asl;
  }
 
  void attach_extra_send_list_function(void (*func)(std::vector<dof_id_type> &, void *),
                                       void * context = nullptr)
  { _extra_send_list_function = func; _extra_send_list_context = context; }
 
  void prepare_send_list ();
 
  void clear_send_list ()
  {
    _send_list.clear();
  }
 
  void reinit_send_list (MeshBase & mesh);
 
 
  const std::vector<dof_id_type> & get_send_list() const { return _send_list; }
 
  const std::vector<dof_id_type> & get_n_nz() const
  {
    libmesh_assert(_sp);
    return _sp->get_n_nz();
  }
 
  const std::vector<dof_id_type> & get_n_oz() const
  {
    libmesh_assert(_sp);
    return _sp->get_n_oz();
  }
 
 
  const SparsityPattern::Build * get_sparsity_pattern() const
  {
    return _sp.get();
  }
 
  unsigned int n_vars() const;
 
  const std::string & variable_name(const unsigned int i) const;
 
  unsigned int n_components(const MeshBase & mesh) const;
 
  bool identify_variable_groups () const;
 
  void identify_variable_groups (const bool);
 
  unsigned int variable_scalar_number (unsigned int var_num,
                                       unsigned int component) const;
 
  const FEType & variable_type (const unsigned int i) const;
 
  const FEType & variable_type (std::string_view var) const;
 
  unsigned int variable_number (std::string_view var) const;
 
  bool has_variable(std::string_view var) const;
 
  void get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const;
 
  unsigned int add_variable (System & sys,
                             std::string_view var,
                             const FEType & type,
                             const std::set<subdomain_id_type> * const active_subdomains = nullptr);
 
  unsigned int add_variables (System & sys,
                              const std::vector<std::string> & vars,
                              const FEType & type,
                              const std::set<subdomain_id_type> * const active_subdomains = nullptr);
 
  unsigned int add_variable_array (System & sys,
                                   const std::vector<std::string> & vars,
                                   const FEType & type,
                                   const std::set<subdomain_id_type> * const active_subdomains = nullptr);
 
  void set_error_on_cyclic_constraint(bool error_on_cyclic_constraint);
  void set_error_on_constraint_loop(bool error_on_constraint_loop);
 
  const VariableGroup & variable_group (const unsigned int c) const;
 
  const Variable & variable (const unsigned int c) const override;
 
  Order variable_order (const unsigned int c) const;
 
  Order variable_group_order (const unsigned int vg) const;
 
  const FEType & variable_group_type (const unsigned int vg) const;
 
  unsigned int n_variable_groups() const
  { return cast_int<unsigned int>(_variable_groups.size()); }
 
  unsigned int n_variables() const override
  { return cast_int<unsigned int>(_variables.size()); }
 
  unsigned int var_group_from_var_number(unsigned int var_num) const;
 
  bool has_blocked_representation() const
  {
    return ((this->n_variable_groups() == 1) && (this->n_variables() > 1));
  }
 
  unsigned int block_size() const
  {
    return (this->has_blocked_representation() ? this->n_variables() : 1);
  }
 
  using DofMapBase::n_dofs;
  dof_id_type n_dofs(const unsigned int vn) const
  {
    dof_id_type n = this->n_local_dofs(vn);
    this->comm().sum(n);
    return n;
  }
 
  dof_id_type n_SCALAR_dofs() const { return _n_SCALAR_dofs; }
 
  using DofMapBase::n_local_dofs;
  dof_id_type n_local_dofs(const unsigned int vn) const
  {
    dof_id_type n;
    this->local_variable_indices(n, _mesh, vn);
    return n;
  }
 
  std::vector<dof_id_type> n_dofs_per_processor(const unsigned int vn) const
  {
    std::vector<dof_id_type> n_local_dofs(this->n_processors(), 0);
    this->comm().allgather(this->n_local_dofs(vn), n_local_dofs);
    return n_local_dofs;
  }
 
  processor_id_type dof_owner(const dof_id_type dof) const
  { std::vector<dof_id_type>::const_iterator ub =
      std::upper_bound(_end_df.begin(), _end_df.end(), dof);
    libmesh_assert (ub != _end_df.end());
    return cast_int<processor_id_type>(ub - _end_df.begin());
  }
 
  void dof_indices (const Elem * const elem,
                    std::vector<dof_id_type> & di) const;
 
  void dof_indices (const Elem * const elem,
                    std::vector<dof_id_type> & di,
                    const unsigned int vn,
                    int p_level = -12345) const override;
 
  void array_dof_indices(const Elem * const elem,
                         std::vector<dof_id_type> & di,
                         const unsigned int vn,
                         int p_level = -12345) const;
 
  void array_dof_indices(const Node * const node,
                         std::vector<dof_id_type> & di,
                         const unsigned int vn) const;
 
  template <typename DofIndicesFunctor>
  void array_dof_indices(const DofIndicesFunctor & functor,
                         std::vector<dof_id_type> & di,
                         const unsigned int vn) const;
 
  template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
  void dof_indices(const Elem * const elem,
                   std::vector<dof_id_type> & di,
                   const unsigned int vn,
                   ScalarDofsFunctor scalar_dofs_functor,
                   FieldDofsFunctor field_dofs_functor,
                   int p_level = -12345) const;
 
  void dof_indices (const Node * const node,
                    std::vector<dof_id_type> & di) const;
 
  void dof_indices (const Node * const node,
                    std::vector<dof_id_type> & di,
                    const unsigned int vn) const override;
 
  void dof_indices (const Elem & elem,
                    unsigned int n,
                    std::vector<dof_id_type> & di,
                    const unsigned int vn) const;
 
#ifdef LIBMESH_ENABLE_AMR
 
  void old_dof_indices (const Elem & elem,
                        unsigned int n,
                        std::vector<dof_id_type> & di,
                        const unsigned int vn) const;
 
#endif // LIBMESH_ENABLE_AMR
 
  void SCALAR_dof_indices (std::vector<dof_id_type> & di,
                           const unsigned int vn,
                           const bool old_dofs=false) const;
 
  bool semilocal_index (dof_id_type dof_index) const;
 
  bool all_semilocal_indices (const std::vector<dof_id_type> & dof_indices) const;
 
  bool local_index (dof_id_type dof_index) const
  { return (dof_index >= this->first_dof()) && (dof_index < this->end_dof()); }
 
  template <typename DofObjectSubclass>
  bool is_evaluable(const DofObjectSubclass & obj,
                    unsigned int var_num = libMesh::invalid_uint) const;
 
  void set_implicit_neighbor_dofs(bool implicit_neighbor_dofs);
 
  void set_verify_dirichlet_bc_consistency(bool val);
 
  bool use_coupled_neighbor_dofs(const MeshBase & mesh) const;
 
  void extract_local_vector (const NumericVector<Number> & Ug,
                             const std::vector<dof_id_type> & dof_indices,
                             DenseVectorBase<Number> & Ue) const;
 
  template <typename T, std::enable_if_t<std::is_same_v<T, dof_id_type> ||
                                         std::is_same_v<T, std::vector<dof_id_type>>, int> = 0>
  void local_variable_indices(T & idx,
                              const MeshBase & mesh,
                              unsigned int var_num) const;
 
  template <typename T,
            std::enable_if_t<std::is_same_v<T, dof_id_type> ||
                                 std::is_same_v<T, std::vector<dof_id_type>>,
                             int> = 0>
  void local_variable_indices(T & idx, unsigned int var_num) const
  { this->local_variable_indices(idx, this->_mesh, var_num); }
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
 
  //--------------------------------------------------------------------
  // Constraint-specific methods
  dof_id_type n_constrained_dofs() const;
 
  dof_id_type n_local_constrained_dofs() const;
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  dof_id_type n_constrained_nodes() const
  { return cast_int<dof_id_type>(_node_constraints.size()); }
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
  void create_dof_constraints (const MeshBase &, Real time=0);
 
  void allgather_recursive_constraints (MeshBase &);
 
  void scatter_constraints (MeshBase &);
 
  void gather_constraints (MeshBase & mesh,
                           std::set<dof_id_type> & unexpanded_dofs,
                           bool look_for_constrainees);
 
  void process_constraints (MeshBase &);
 
  void check_for_cyclic_constraints();
  void check_for_constraint_loops();
 
  void add_constraint_row (const dof_id_type dof_number,
                           const DofConstraintRow & constraint_row,
                           const Number constraint_rhs,
                           const bool forbid_constraint_overwrite);
 
  void add_adjoint_constraint_row (const unsigned int qoi_index,
                                   const dof_id_type dof_number,
                                   const DofConstraintRow & constraint_row,
                                   const Number constraint_rhs,
                                   const bool forbid_constraint_overwrite);
 
  void add_constraint_row (const dof_id_type dof_number,
                           const DofConstraintRow & constraint_row,
                           const bool forbid_constraint_overwrite = true)
  { add_constraint_row(dof_number, constraint_row, 0., forbid_constraint_overwrite); }
 
  DofConstraints::const_iterator constraint_rows_begin() const
  { return _dof_constraints.begin(); }
 
  DofConstraints::const_iterator constraint_rows_end() const
  { return _dof_constraints.end(); }
 
  const DofConstraints & get_dof_constraints() const { return _dof_constraints; }
 
  void stash_dof_constraints()
  {
    libmesh_assert(_stashed_dof_constraints.empty());
    _dof_constraints.swap(_stashed_dof_constraints);
  }
 
  void unstash_dof_constraints()
  {
    libmesh_assert(_dof_constraints.empty());
    _dof_constraints.swap(_stashed_dof_constraints);
  }
 
  void swap_dof_constraints()
  {
    _dof_constraints.swap(_stashed_dof_constraints);
  }
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  NodeConstraints::const_iterator node_constraint_rows_begin() const
  { return _node_constraints.begin(); }
 
  NodeConstraints::const_iterator node_constraint_rows_end() const
  { return _node_constraints.end(); }
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
  bool is_constrained_dof (const dof_id_type dof) const;
 
  bool has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const;
 
  bool has_heterogenous_adjoint_constraints (const unsigned int qoi_num) const
  {
    return this->has_heterogeneous_adjoint_constraints (qoi_num);
  }
 
  Number has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
                                               const dof_id_type dof) const;
 
  Number has_heterogenous_adjoint_constraint (const unsigned int qoi_num,
                                              const dof_id_type dof) const
  {
    return this->has_heterogeneous_adjoint_constraint (qoi_num, dof);
  }
 
  DofConstraintValueMap & get_primal_constraint_values();
 
  bool is_constrained_node (const Node * node) const;
 
  void print_dof_constraints(std::ostream & os=libMesh::out,
                             bool print_nonlocal=false) const;
 
  std::string get_local_constraints(bool print_nonlocal=false) const;
 
 
  std::pair<Real, Real> max_constraint_error(const System & system,
                                             NumericVector<Number> * v = nullptr) const;
 
#endif // LIBMESH_ENABLE_CONSTRAINTS
 
  //--------------------------------------------------------------------
  // Constraint-specific methods
  // Some of these methods are enabled (but inlined away to nothing)
  // when constraints are disabled at configure-time.  This is to
  // increase API compatibility of user code with different library
  // builds.
 
  void constrain_element_matrix (DenseMatrix<Number> & matrix,
                                 std::vector<dof_id_type> & elem_dofs,
                                 bool asymmetric_constraint_rows = true) const;
 
  void constrain_element_matrix (DenseMatrix<Number> & matrix,
                                 std::vector<dof_id_type> & row_dofs,
                                 std::vector<dof_id_type> & col_dofs,
                                 bool asymmetric_constraint_rows = true) const;
 
  void constrain_element_vector (DenseVector<Number> & rhs,
                                 std::vector<dof_id_type> & dofs,
                                 bool asymmetric_constraint_rows = true) const;
 
  void constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
                                            DenseVector<Number> & rhs,
                                            std::vector<dof_id_type> & elem_dofs,
                                            bool asymmetric_constraint_rows = true) const;
 
  void heterogeneously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
                                                            DenseVector<Number> & rhs,
                                                            std::vector<dof_id_type> & elem_dofs,
                                                            bool asymmetric_constraint_rows = true,
                                                            int qoi_index = -1) const;
 
  /*
   * Backwards compatibility with misspelling.
   */
  void heterogenously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
                                                           DenseVector<Number> & rhs,
                                                           std::vector<dof_id_type> & elem_dofs,
                                                           bool asymmetric_constraint_rows = true,
                                                           int qoi_index = -1) const
  {
    return this->heterogeneously_constrain_element_matrix_and_vector
      (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
  }
 
  void heterogeneously_constrain_element_vector (const DenseMatrix<Number> & matrix,
                                                 DenseVector<Number> & rhs,
                                                 std::vector<dof_id_type> & elem_dofs,
                                                 bool asymmetric_constraint_rows = true,
                                                 int qoi_index = -1) const;
 
  /*
   * Backwards compatibility with misspelling.
   */
  void heterogenously_constrain_element_vector (const DenseMatrix<Number> & matrix,
                                                DenseVector<Number> & rhs,
                                                std::vector<dof_id_type> & elem_dofs,
                                                bool asymmetric_constraint_rows = true,
                                                int qoi_index = -1) const
  {
    return this->heterogeneously_constrain_element_vector
      (matrix, rhs, elem_dofs, asymmetric_constraint_rows, qoi_index);
  }
 
  void heterogeneously_constrain_element_jacobian_and_residual (DenseMatrix<Number> & matrix,
                                                                DenseVector<Number> & rhs,
                                                                std::vector<dof_id_type> & elem_dofs,
                                                                NumericVector<Number> & solution_local) const;
 
  void heterogeneously_constrain_element_residual (DenseVector<Number> & rhs,
                                                   std::vector<dof_id_type> & elem_dofs,
                                                   NumericVector<Number> & solution_local) const;
 
 
  void constrain_element_residual (DenseVector<Number> & rhs,
                                   std::vector<dof_id_type> & elem_dofs,
                                   NumericVector<Number> & solution_local) const;
 
  void constrain_element_dyad_matrix (DenseVector<Number> & v,
                                      DenseVector<Number> & w,
                                      std::vector<dof_id_type> & row_dofs,
                                      bool asymmetric_constraint_rows = true) const;
 
  void constrain_nothing (std::vector<dof_id_type> & dofs) const;
 
  void enforce_constraints_exactly (const System & system,
                                    NumericVector<Number> * v = nullptr,
                                    bool homogeneous = false) const;
 
  void enforce_adjoint_constraints_exactly (NumericVector<Number> & v,
                                            unsigned int q) const;
 
  void enforce_constraints_on_residual (const NonlinearImplicitSystem & system,
                                        NumericVector<Number> * rhs,
                                        NumericVector<Number> const * solution,
                                        bool homogeneous = true) const;
 
  void enforce_constraints_on_jacobian (const NonlinearImplicitSystem & system,
                                        SparseMatrix<Number> * jac) const;
 
#ifdef LIBMESH_ENABLE_PERIODIC
 
  //--------------------------------------------------------------------
  // PeriodicBoundary-specific methods
 
  void add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary);
 
  void add_periodic_boundary (const PeriodicBoundaryBase & boundary, const PeriodicBoundaryBase & inverse_boundary);
 
  bool is_periodic_boundary (const boundary_id_type boundaryid) const;
 
  PeriodicBoundaries * get_periodic_boundaries()
  {
    return _periodic_boundaries.get();
  }
 
  const PeriodicBoundaries * get_periodic_boundaries() const
  {
    return _periodic_boundaries.get();
  }
 
#endif // LIBMESH_ENABLE_PERIODIC
 
 
#ifdef LIBMESH_ENABLE_DIRICHLET
 
  //--------------------------------------------------------------------
  // DirichletBoundary-specific methods
 
  void add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
 
  void add_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
                                       unsigned int q);
 
  void remove_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary);
 
  void remove_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
                                          unsigned int q);
 
  const DirichletBoundaries * get_dirichlet_boundaries() const
  {
    return _dirichlet_boundaries.get();
  }
 
  DirichletBoundaries * get_dirichlet_boundaries()
  {
    return _dirichlet_boundaries.get();
  }
 
  bool has_adjoint_dirichlet_boundaries(unsigned int q) const;
 
  const DirichletBoundaries *
  get_adjoint_dirichlet_boundaries(unsigned int q) const;
 
  DirichletBoundaries *
  get_adjoint_dirichlet_boundaries(unsigned int q);
 
  void check_dirichlet_bcid_consistency (const MeshBase & mesh,
                                         const DirichletBoundary & boundary) const;
#endif // LIBMESH_ENABLE_DIRICHLET
 
 
#ifdef LIBMESH_ENABLE_AMR
 
  //--------------------------------------------------------------------
  // AMR-specific methods
 
  // void augment_send_list_for_projection(const MeshBase &);
 
#ifdef LIBMESH_ENABLE_AMR
 
  void old_dof_indices (const Elem * const elem,
                        std::vector<dof_id_type> & di,
                        const unsigned int vn = libMesh::invalid_uint) const;
 
#endif // LIBMESH_ENABLE_AMR
 
  void constrain_p_dofs (unsigned int var,
                         const Elem * elem,
                         unsigned int s,
                         unsigned int p);
 
#endif // LIBMESH_ENABLE_AMR
 
  void reinit
    (MeshBase & mesh,
     const std::map<const Node *, std::set<subdomain_id_type>> &
       constraining_subdomains);
 
  virtual void clear () override;
 
  void print_info(std::ostream & os=libMesh::out) const;
 
  std::string get_info() const;
 
  CouplingMatrix * _dof_coupling;
 
  unsigned int sys_number() const;
 
  std::unique_ptr<SparsityPattern::Build> build_sparsity(const MeshBase & mesh,
                                                         bool calculate_constrained = false,
                                                         bool use_condensed_system = false) const;
 
  void should_p_refine(unsigned int g, bool p_refine);
 
  bool should_p_refine(unsigned int g) const;
 
  bool should_p_refine_var(unsigned int var) const;
 
  // Prevent bad user implicit conversions
  void should_p_refine(FEFamily, bool) = delete;
  void should_p_refine(Order, bool) = delete;
  bool should_p_refine(FEFamily) const = delete;
  bool should_p_refine(Order) const = delete;
 
  void create_static_condensation(MeshBase & mesh, System & system);
 
  bool has_static_condensation() const { return _sc.get(); }
 
  StaticCondensationDofMap & get_static_condensation();
 
  const StaticCondensationDofMap & get_static_condensation() const;
 
  void reinit_static_condensation();
 
private:
 
  const std::pair<unsigned int, unsigned int> &
  get_variable_array(unsigned int vi) const;
 
  void _dof_indices (const Elem & elem,
                     int p_level,
                     std::vector<dof_id_type> & di,
                     const unsigned int vg,
                     const unsigned int vig,
                     const Node * const * nodes,
                     unsigned int       n_nodes,
                     const unsigned int v
#ifdef DEBUG
                     ,
                     std::size_t & tot_size
#endif
                    ) const;
 
  template <typename FieldDofsFunctor>
  void _dof_indices (const Elem & elem,
                     int p_level,
                     std::vector<dof_id_type> & di,
                     const unsigned int vg,
                     const unsigned int vig,
                     const Node * const * nodes,
                     unsigned int       n_nodes,
                     const unsigned int v,
#ifdef DEBUG
                     std::size_t & tot_size,
#endif
                     FieldDofsFunctor field_dofs_functor) const;
 
  void _node_dof_indices (const Elem & elem,
                          unsigned int n,
                          const DofObject & obj,
                          std::vector<dof_id_type> & di,
                          const unsigned int vn) const;
 
  void invalidate_dofs(MeshBase & mesh) const;
 
  DofObject * node_ptr(MeshBase & mesh, dof_id_type i) const;
 
  DofObject * elem_ptr(MeshBase & mesh, dof_id_type i) const;
 
  typedef DofObject * (DofMap::*dofobject_accessor)
    (MeshBase & mesh, dof_id_type i) const;
 
  template<typename iterator_type>
  void set_nonlocal_dof_objects(iterator_type objects_begin,
                                iterator_type objects_end,
                                MeshBase & mesh,
                                dofobject_accessor objects);
 
  std::map<const Node *, std::set<subdomain_id_type>>
  calculate_constraining_subdomains();
 
  void distribute_local_dofs_var_major
    (dof_id_type & next_free_dof,
     MeshBase & mesh,
     const std::map<const Node *, std::set<subdomain_id_type>> &
       constraining_subdomains);
 
  void distribute_local_dofs_node_major
    (dof_id_type & next_free_dof,
     MeshBase & mesh,
     const std::map<const Node *, std::set<subdomain_id_type>> &
       constraining_subdomains);
 
  /*
   * Helper method for the above two to count + distriubte SCALAR dofs
   */
  void distribute_scalar_dofs (dof_id_type & next_free_dof);
 
#ifdef DEBUG
  /*
   * Internal assertions for distribute_local_dofs_*
   */
  void assert_no_nodes_missed(MeshBase & mesh);
#endif
 
  /*
   * A utility method for obtaining a set of elements to ghost along
   * with merged coupling matrices.
   */
  typedef std::set<std::unique_ptr<CouplingMatrix>, Utility::CompareUnderlying> CouplingMatricesSet;
  static void
  merge_ghost_functor_outputs (GhostingFunctor::map_type & elements_to_ghost,
                               CouplingMatricesSet & temporary_coupling_matrices,
                               const GhostingFunctorIterator & gf_begin,
                               const GhostingFunctorIterator & gf_end,
                               const MeshBase::const_element_iterator & elems_begin,
                               const MeshBase::const_element_iterator & elems_end,
                               processor_id_type p);
 
  void add_neighbors_to_send_list(MeshBase & mesh);
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
 
  void build_constraint_matrix (DenseMatrix<Number> & C,
                                std::vector<dof_id_type> & elem_dofs,
                                const bool called_recursively=false) const;
 
  void build_constraint_matrix_and_vector (DenseMatrix<Number> & C,
                                           DenseVector<Number> & H,
                                           std::vector<dof_id_type> & elem_dofs,
                                           int qoi_index = -1,
                                           const bool called_recursively=false) const;
 
  void find_connected_dofs (std::vector<dof_id_type> & elem_dofs) const;
 
  void find_connected_dof_objects (std::vector<const DofObject *> & objs) const;
 
  void add_constraints_to_send_list();
 
  void process_mesh_constraint_rows(const MeshBase & mesh);
 
#endif // LIBMESH_ENABLE_CONSTRAINTS
 
  bool _error_on_constraint_loop;
 
  bool _constrained_sparsity_construction;
 
  std::vector<Variable> _variables;
 
  std::vector<VariableGroup> _variable_groups;
 
  std::vector<unsigned int> _variable_group_numbers;
 
  std::unordered_map<unsigned int, unsigned int> _var_to_vg;
 
  std::map<std::string, unsigned int, std::less<>> _variable_numbers;
 
  std::vector<std::pair<unsigned int, unsigned int>> _array_variables;
 
  bool _identify_variable_groups = true;
 
  const unsigned int _sys_number;
 
  MeshBase & _mesh;
 
  std::vector<SparseMatrix<Number> * > _matrices;
 
  std::vector<dof_id_type> _first_scalar_df;
 
  std::vector<dof_id_type> _send_list;
 
  SparsityPattern::AugmentSparsityPattern * _augment_sparsity_pattern;
 
  void (*_extra_sparsity_function)(SparsityPattern::Graph &,
                                   std::vector<dof_id_type> & n_nz,
                                   std::vector<dof_id_type> & n_oz,
                                   void *);
  void * _extra_sparsity_context;
 
  AugmentSendList * _augment_send_list;
 
  void (*_extra_send_list_function)(std::vector<dof_id_type> &, void *);
 
  void * _extra_send_list_context;
 
  std::unique_ptr<DefaultCoupling> _default_coupling;
 
  std::unique_ptr<DefaultCoupling> _default_evaluating;
 
  std::vector<GhostingFunctor *> _algebraic_ghosting_functors;
 
  std::vector<GhostingFunctor *> _coupling_functors;
 
  std::map<GhostingFunctor *, std::shared_ptr<GhostingFunctor> > _shared_functors;
 
  bool need_full_sparsity_pattern;
 
  std::unique_ptr<SparsityPattern::Build> _sp;
 
  dof_id_type _n_SCALAR_dofs;
 
#ifdef LIBMESH_ENABLE_AMR
 
  std::vector<dof_id_type> _first_old_scalar_df;
#endif
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
  DofConstraints _dof_constraints, _stashed_dof_constraints;
 
  DofConstraintValueMap      _primal_constraint_values;
 
  AdjointDofConstraintValues _adjoint_constraint_values;
#endif
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  NodeConstraints _node_constraints;
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
 
#ifdef LIBMESH_ENABLE_PERIODIC
  std::unique_ptr<PeriodicBoundaries> _periodic_boundaries;
#endif
 
#ifdef LIBMESH_ENABLE_DIRICHLET
  std::unique_ptr<DirichletBoundaries> _dirichlet_boundaries;
 
  std::vector<std::unique_ptr<DirichletBoundaries>> _adjoint_dirichlet_boundaries;
#endif
 
  friend class SparsityPattern::Build;
 
  bool _implicit_neighbor_dofs_initialized;
  bool _implicit_neighbor_dofs;
 
  bool _verify_dirichlet_bc_consistency;
 
  std::unique_ptr<StaticCondensationDofMap> _sc;
};
 
 
// ------------------------------------------------------------
// Dof Map inline member functions
inline
unsigned int DofMap::sys_number() const
{
  return _sys_number;
}
 
 
 
inline
const VariableGroup & DofMap::variable_group (const unsigned int g) const
{
  libmesh_assert_less (g, _variable_groups.size());
 
  return _variable_groups[g];
}
 
 
 
inline
const Variable & DofMap::variable (const unsigned int c) const
{
  libmesh_assert_less (c, _variables.size());
 
  return _variables[c];
}
 
 
 
inline
Order DofMap::variable_order (const unsigned int c) const
{
  libmesh_assert_less (c, _variables.size());
 
  return _variables[c].type().order;
}
 
 
 
inline
Order DofMap::variable_group_order (const unsigned int vg) const
{
  libmesh_assert_less (vg, _variable_groups.size());
 
  return _variable_groups[vg].type().order;
}
 
 
 
inline
const FEType & DofMap::variable_type (const unsigned int c) const
{
  libmesh_assert_less (c, _variables.size());
 
  return _variables[c].type();
}
 
 
 
inline
const FEType & DofMap::variable_group_type (const unsigned int vg) const
{
  libmesh_assert_less (vg, _variable_groups.size());
 
  return _variable_groups[vg].type();
}
 
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
 
 
inline
bool DofMap::is_constrained_node (const Node *
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
                                  node
#endif
                                  ) const
{
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  if (_node_constraints.count(node))
    return true;
#endif
 
  return false;
}
 
 
inline
bool DofMap::is_constrained_dof (const dof_id_type dof) const
{
  if (_dof_constraints.count(dof))
    return true;
 
  return false;
}
 
 
inline
bool DofMap::has_heterogeneous_adjoint_constraints (const unsigned int qoi_num) const
{
  AdjointDofConstraintValues::const_iterator it =
    _adjoint_constraint_values.find(qoi_num);
  if (it == _adjoint_constraint_values.end())
    return false;
  if (it->second.empty())
    return false;
 
  return true;
}
 
 
inline
Number DofMap::has_heterogeneous_adjoint_constraint (const unsigned int qoi_num,
                                                    const dof_id_type dof) const
{
  AdjointDofConstraintValues::const_iterator it =
    _adjoint_constraint_values.find(qoi_num);
  if (it != _adjoint_constraint_values.end())
    {
      DofConstraintValueMap::const_iterator rhsit =
        it->second.find(dof);
      if (rhsit == it->second.end())
        return 0;
      else
        return rhsit->second;
    }
 
  return 0;
}
 
 
 
inline
DofConstraintValueMap & DofMap::get_primal_constraint_values()
{
  return _primal_constraint_values;
}
 
 
 
#else
 
//--------------------------------------------------------------------
// Constraint-specific methods get inlined into nothing if
// constraints are disabled, so there's no reason for users not to
// use them.
 
inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
                                              std::vector<dof_id_type> &,
                                              bool) const {}
 
inline void DofMap::constrain_element_matrix (DenseMatrix<Number> &,
                                              std::vector<dof_id_type> &,
                                              std::vector<dof_id_type> &,
                                              bool) const {}
 
inline void DofMap::constrain_element_vector (DenseVector<Number> &,
                                              std::vector<dof_id_type> &,
                                              bool) const {}
 
inline void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number> &,
                                                         DenseVector<Number> &,
                                                         std::vector<dof_id_type> &,
                                                         bool) const {}
 
inline void DofMap::heterogeneously_constrain_element_matrix_and_vector
  (DenseMatrix<Number> &, DenseVector<Number> &,
   std::vector<dof_id_type> &, bool, int) const {}
 
inline void DofMap::heterogeneously_constrain_element_vector
  (const DenseMatrix<Number> &, DenseVector<Number> &,
   std::vector<dof_id_type> &, bool, int) const {}
 
inline void DofMap::constrain_element_dyad_matrix (DenseVector<Number> &,
                                                   DenseVector<Number> &,
                                                   std::vector<dof_id_type> &,
                                                   bool) const {}
 
inline void DofMap::constrain_nothing (std::vector<dof_id_type> &) const {}
 
inline void DofMap::enforce_constraints_exactly (const System &,
                                                 NumericVector<Number> *,
                                                 bool) const {}
 
inline void DofMap::enforce_adjoint_constraints_exactly (NumericVector<Number> &,
                                                         unsigned int) const {}
 
 
inline void DofMap::enforce_constraints_on_residual
  (const NonlinearImplicitSystem &,
   NumericVector<Number> *,
   NumericVector<Number> const *,
   bool) const {}
 
inline void DofMap::enforce_constraints_on_jacobian
  (const NonlinearImplicitSystem &,
   SparseMatrix<Number> *) const {}
 
#endif // LIBMESH_ENABLE_CONSTRAINTS
 
 
 
inline
void DofMap::set_constrained_sparsity_construction(bool use_constraints)
{
  // This got only partly finished...
  if (use_constraints)
    libmesh_not_implemented();
 
#ifdef LIBMESH_ENABLE_CONSTRAINTS
  _constrained_sparsity_construction = use_constraints;
#endif
  libmesh_ignore(use_constraints);
}
 
inline
void DofMap::full_sparsity_pattern_needed()
{
  need_full_sparsity_pattern = true;
}
 
inline
bool DofMap::constrained_sparsity_construction()
{
#ifdef LIBMESH_ENABLE_CONSTRAINTS
  return _constrained_sparsity_construction;
#else
  return true;
#endif
}
 
inline
void DofMap::should_p_refine(const unsigned int g, const bool p_refine)
{
#ifdef LIBMESH_ENABLE_AMR
  VariableGroup & var = _variable_groups[g];
  var.type().p_refinement = p_refine;
 
  for (auto v : make_range(var.first_scalar_number(0),
                           var.first_scalar_number(0) +
                           var.n_variables()))
    this->_variables[v].type().p_refinement = p_refine;
 
 
#else
  libmesh_ignore(g, p_refine);
#endif
}
 
inline
bool DofMap::should_p_refine(const unsigned int g) const
{
#ifdef LIBMESH_ENABLE_AMR
  const VariableGroup & var = this->variable_group(g);
  return var.type().p_refinement;
#else
  libmesh_ignore(g);
  return false;
#endif
}
 
inline
unsigned int DofMap::var_group_from_var_number(const unsigned int var_num) const
{
  libmesh_assert(var_num < n_variables());
  return libmesh_map_find(_var_to_vg, var_num);
}
 
inline
bool DofMap::should_p_refine_var(const unsigned int var) const
{
#ifdef LIBMESH_ENABLE_AMR
  const auto vg = this->var_group_from_var_number(var);
  return this->should_p_refine(vg);
#else
  libmesh_ignore(var);
  return false;
#endif
}
 
template <typename FieldDofsFunctor>
void DofMap::_dof_indices (const Elem & elem,
                           int p_level,
                           std::vector<dof_id_type> & di,
                           const unsigned int vg,
                           const unsigned int vig,
                           const Node * const * nodes,
                           unsigned int       n_nodes,
                           const unsigned int v,
#ifdef DEBUG
                           std::size_t & tot_size,
#endif
                           FieldDofsFunctor field_dofs_functor) const
{
  const VariableGroup & var = this->variable_group(vg);
 
  if (var.active_on_subdomain(elem.subdomain_id()))
    {
      const ElemType type        = elem.type();
      const unsigned int sys_num = this->sys_number();
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
      const bool is_inf          = elem.infinite();
#endif
 
      const bool extra_hanging_dofs =
        FEInterface::extra_hanging_dofs(var.type());
 
      FEType fe_type = var.type();
 
      const bool add_p_level = fe_type.p_refinement;
 
#ifdef DEBUG
      // The number of dofs per element is non-static for subdivision FE
      if (var.type().family == SUBDIVISION)
        tot_size += n_nodes;
      else
        // FIXME: Is the passed-in p_level just elem.p_level()? If so,
        // this seems redundant.
        tot_size += FEInterface::n_dofs(fe_type, add_p_level*p_level, &elem);
#endif
 
      // The total Order is not required when getting the function
      // pointer, it is only needed when the function is called (see
      // below).
      const FEInterface::n_dofs_at_node_ptr ndan =
        FEInterface::n_dofs_at_node_function(fe_type, &elem);
 
      // Get the node-based DOF numbers
      for (unsigned int n=0; n != n_nodes; n++)
        {
          const Node & node = *nodes[n];
 
          // Cache the intermediate lookups that are common to every
          // component
#ifdef DEBUG
          const std::pair<unsigned int, unsigned int>
            vg_and_offset = node.var_to_vg_and_offset(sys_num,v);
          libmesh_assert_equal_to (vg, vg_and_offset.first);
          libmesh_assert_equal_to (vig, vg_and_offset.second);
#endif
          const unsigned int n_comp = node.n_comp_group(sys_num,vg);
 
          // There is a potential problem with h refinement.  Imagine a
          // quad9 that has a linear FE on it.  Then, on the hanging side,
          // it can falsely identify a DOF at the mid-edge node. This is why
          // we go through FEInterface instead of node.n_comp() directly.
          const unsigned int nc =
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
            is_inf ?
            FEInterface::n_dofs_at_node(fe_type, add_p_level*p_level, &elem, n) :
#endif
            ndan (type, fe_type.order + add_p_level*p_level, n);
 
          // If this is a non-vertex on a hanging node with extra
          // degrees of freedom, we use the non-vertex dofs (which
          // come in reverse order starting from the end, to
          // simplify p refinement)
          if (extra_hanging_dofs && !elem.is_vertex(n))
            {
              const int dof_offset = n_comp - nc;
 
              // We should never have fewer dofs than necessary on a
              // node unless we're getting indices on a parent element,
              // and we should never need the indices on such a node
              if (dof_offset < 0)
                {
                  libmesh_assert(!elem.active());
                  di.resize(di.size() + nc, DofObject::invalid_id);
                }
              else
                for (int i=int(n_comp)-1; i>=dof_offset; i--)
                  {
                    const dof_id_type d =
                      node.dof_number(sys_num, vg, vig, i, n_comp);
                    libmesh_assert_not_equal_to (d, DofObject::invalid_id);
                    field_dofs_functor(elem, n, v, di, d);
                  }
            }
          // If this is a vertex or an element without extra hanging
          // dofs, our dofs come in forward order coming from the
          // beginning
          else
            {
              // We have a good component index only if it's being
              // used on this FE type (nc) *and* it's available on
              // this DofObject (n_comp).
              const unsigned int good_nc = std::min(n_comp, nc);
              for (unsigned int i=0; i!=good_nc; ++i)
                {
                  const dof_id_type d =
                    node.dof_number(sys_num, vg, vig, i, n_comp);
                  libmesh_assert_not_equal_to (d, DofObject::invalid_id);
                  libmesh_assert_less (d, this->n_dofs());
                  field_dofs_functor(elem, n, v, di, d);
                }
 
              // With fewer good component indices than we need, e.g.
              // due to subdomain expansion, the remaining expected
              // indices are marked invalid.
              if (n_comp < nc)
                for (unsigned int i=n_comp; i!=nc; ++i)
                  di.push_back(DofObject::invalid_id);
            }
        }
 
      // If there are any element-based DOF numbers, get them
      const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, add_p_level*p_level, &elem);
 
      // We should never have fewer dofs than necessary on an
      // element unless we're getting indices on a parent element
      // (and we should never need those indices) or off-domain for a
      // subdomain-restricted variable (where invalid_id is the
      // correct thing to return)
      if (nc != 0)
        {
          const unsigned int n_comp = elem.n_comp_group(sys_num,vg);
          if (elem.n_systems() > sys_num && nc <= n_comp)
            {
              for (unsigned int i=0; i<nc; i++)
                {
                  const dof_id_type d =
                    elem.dof_number(sys_num, vg, vig, i, n_comp);
                  libmesh_assert_not_equal_to (d, DofObject::invalid_id);
 
                  field_dofs_functor(elem, invalid_uint, v, di, d);
                }
            }
          else
            {
              libmesh_assert(!elem.active() || fe_type.family == LAGRANGE || fe_type.family == SUBDIVISION);
              di.resize(di.size() + nc, DofObject::invalid_id);
            }
        }
    }
}
 
 
 
template <typename ScalarDofsFunctor, typename FieldDofsFunctor>
void DofMap::dof_indices (const Elem * const elem,
                          std::vector<dof_id_type> & di,
                          const unsigned int vn,
                          ScalarDofsFunctor scalar_dofs_functor,
                          FieldDofsFunctor field_dofs_functor,
                          int p_level) const
{
  // We now allow elem==nullptr to request just SCALAR dofs
  // libmesh_assert(elem);
 
  // dof_indices() is a relatively light-weight function that is
  // called millions of times in normal codes. Therefore, it is not a
  // good candidate for logging, since the cost of the logging code
  // itself is roughly on par with the time required to call
  // dof_indices().
  // LOG_SCOPE("dof_indices()", "DofMap");
 
  // Clear the DOF indices vector
  di.clear();
 
  // Use the default p refinement level?
  if (p_level == -12345)
    p_level = elem ? elem->p_level() : 0;
 
  const unsigned int vg = this->_variable_group_numbers[vn];
  const VariableGroup & var = this->variable_group(vg);
  const unsigned int vig = vn - var.number();
 
#ifdef DEBUG
  // Check that sizes match in DEBUG mode
  std::size_t tot_size = 0;
#endif
 
  if (elem && elem->type() == TRI3SUBDIVISION)
    {
      // Subdivision surface FE require the 1-ring around elem
      const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
 
      // Ghost subdivision elements have no real dofs
      if (!sd_elem->is_ghost())
        {
          // Determine the nodes contributing to element elem
          std::vector<const Node *> elem_nodes;
          MeshTools::Subdivision::find_one_ring(sd_elem, elem_nodes);
 
          _dof_indices(*elem, p_level, di, vg, vig, elem_nodes.data(),
                       cast_int<unsigned int>(elem_nodes.size()), vn,
#ifdef DEBUG
                       tot_size,
#endif
                       field_dofs_functor);
        }
 
      return;
    }
 
  // Get the dof numbers
  if (var.type().family == SCALAR &&
      (!elem ||
       var.active_on_subdomain(elem->subdomain_id())))
    {
#ifdef DEBUG
      tot_size += var.type().order;
#endif
      std::vector<dof_id_type> di_new;
      this->SCALAR_dof_indices(di_new,vn);
      scalar_dofs_functor(*elem, di, di_new);
    }
  else if (elem)
    _dof_indices(*elem, p_level, di, vg, vig, elem->get_nodes(),
                 elem->n_nodes(), vn,
#ifdef DEBUG
                 tot_size,
#endif
                 field_dofs_functor);
 
#ifdef DEBUG
  libmesh_assert_equal_to (tot_size, di.size());
#endif
}
 
inline
StaticCondensationDofMap & DofMap::get_static_condensation()
{
  libmesh_assert(_sc);
  return *_sc;
}
 
inline
const StaticCondensationDofMap & DofMap::get_static_condensation() const
{
  libmesh_assert(_sc);
  return *_sc;
}
 
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
 
#endif // LIBMESH_DOF_MAP_H