* K C \f$ and \f$ C^T F \f$ on their unconstrained rows and columns, so dropping their constrained
 * rows and columns at insertion assembles the reduced system exactly.
 */
class ConstraintEliminationDofMap : public ReducedDofMap
{
public:
  ConstraintEliminationDofMap(MeshBase & mesh, System & system, const DofMap & dof_map);

  /**
   * Checks whether we eliminate constrained degrees of freedom
   */
  bool has_constraint_elimination() const { return _constraint_elimination.get(); }

  /**
   * @returns The map of the unconstrained degrees of freedom. This should have been already added

}

inline
const ConstraintEliminationDofMap & DofMap::get_constraint_elimination() const
{
  libmesh_assert(_constraint_elimination);
  return *_constraint_elimination;
}

inline const std::pair<unsigned int, unsigned int> &

  // (*) recording of requested variables happens in '_request_variable()'
  const variable* sv = _request_variable(VarName);
  if (sv == 0)
    return Default;

  const std::string*  element = sv->get_element(Idx);
  if (element == 0)

 * space. This is the information needed to assemble reduced matrices and to build field splits
 * of the reduced space.
 */
class ReducedDofMap : public DofMapBase
{
public:
  virtual ~ReducedDofMap();

  /**
   * Whether \p reinit() has built the reduced space
   */
  bool initialized() const { return _initialized; }

  virtual unsigned int n_variables() const override;

   * \returns The local retained degrees of freedom in the full numbering, sorted. The reduced
   * index of entry \p i is <tt>first_dof() + i</tt>
   */
  const std::vector<dof_id_type> & local_retained_dofs() const { return _local_retained_dofs; }

  /**
   * Nests this reduced space in \p parent's, or in the full space if \p parent is null. Only

   * this one reduces: its vectors are the ones operators on this space are prolonged to and
   * restricted from. Takes effect at the next \p reinit(), which must follow the parent's
   */
  void set_parent(const ReducedDofMap * parent) { _parent = parent; }

  /**
   * \returns The reduced space this one is nested in, or null for the full space
   */
  const ReducedDofMap * parent() const { return _parent; }

  /**
   * \returns The index of \p full_dof in the parent space, or \p DofObject::invalid_id if the

   * \returns The local retained degrees of freedom in the parent numbering, sorted. Entry \p i is
   * the parent index of reduced index <tt>first_dof() + i</tt>
   */
  const std::vector<dof_id_type> & local_retained_parent_dofs() const
  { return _local_retained_parent_dofs; }

  /**
   * \returns The nonlocal reduced degrees of freedom used by active local elements, sorted: the
   * ghost entries of a ghosted vector on this space
   */
  const std::vector<dof_id_type> & send_list() const { return _send_list; }

  /**
   * \returns The ghost entries of a ghosted vector on the parent space

  /**
   * \returns The number of on-processor nonzeros in each local row of the reduced space
   */
  const std::vector<dof_id_type> & reduced_n_nz() const { return _reduced_nnz; }

  /**
   * \returns The number of off-processor nonzeros in each local row of the reduced space
   */
  const std::vector<dof_id_type> & reduced_n_oz() const { return _reduced_noz; }

  /**
   * \returns The reduced index of \p full_dof, or \p DofObject::invalid_id if \p full_dof is

   * This space's rule for retaining degree of freedom \p scalar_dof of a SCALAR variable, among
   * those the parent retains
   */
  virtual bool retain_scalar_dof(dof_id_type /*scalar_dof*/) const { return true; }

protected:
  /**

  return *_reduced_system;
}

inline const SparsityPattern::Build & ReducedDofMap::reduced_sparsity() const
{
  libmesh_assert(_reduced_sp);
  return *_reduced_sp;
}

}

inline
bool Elem::operator != (const Elem & rhs) const
{
  return !(*this == rhs);
}

inline
const Elem * Elem::raw_child_ptr (unsigned int i) const
{
  if (!_children)
    return nullptr;

  return _children[i];
}

inline

  if (elem->has_children())
    for (auto & c : elem->child_ref_range())
      if (!c.is_remote() && c.has_neighbor(neighbor_in))
        ElemInternal::total_family_tree_by_neighbor (&c, family, neighbor_in, false);
}

template <class T>
void
total_family_tree_by_subneighbor(T elem,
                                 std::vector<T> & family,
                                 T neighbor_in,
                                 T subneighbor,
                                 bool reset = true)
{
  // Clear the vector if the flag reset tells us to.
  if (reset)
    family.clear();

  // To simplify this function we need an existing neighbor

  libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));

  // Add this element to the family tree if applicable.
  if (neighbor_in == subneighbor)
    family.push_back(elem);

  // Recurse into the elements children, if it has any.
  // Do not clear the vector any more.
  if (elem->has_children())
    for (auto & c : elem->child_ref_range())
      if (!c.is_remote())
        for (auto child_neigh : c.neighbor_ptr_range())
          if (child_neigh &&
              (child_neigh == neighbor_in ||
               (child_neigh->parent() == neighbor_in && child_neigh->is_ancestor_of(subneighbor))))
            c.total_family_tree_by_subneighbor(family, child_neigh, subneighbor, false);
}

#ifdef LIBMESH_ENABLE_PERIODIC
  // Set PeriodicBoundaries to couple
  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs) override
  { _periodic_bcs = periodic_bcs; }
#endif

  /**

   * may want to delete the no-longer-relevant parts of those caches
   * after a redistribution is complete.
   */
  virtual void delete_remote_elements () {};

protected:
  const MeshBase * _mesh;

  virtual void redistribute () override
  { this->mesh_reinit(); }

  virtual void delete_remote_elements() override
  { this->mesh_reinit(); }

  /**
   * For the specified range of active elements, find the elements

  virtual void redistribute () override
  { this->mesh_reinit(); }

  virtual void delete_remote_elements() override
  { this->mesh_reinit(); }

  /**
   * For the specified range of active elements, find the elements

  /**
   * Default ghosting functor
   */
  GhostingFunctor & default_ghosting() { return *_default_ghosting; }

  /**
   * Constructs a list of all subdomain identifiers in the local mesh if

   * this mesh as a distinct lower-dimensional boundary (or boundary
   * subset) mesh, the original mesh will be returned here.
   */
  const MeshBase & interior_mesh() const { return *_interior_mesh; }

  /**
   * \return A writeable reference to the interior mesh.

    libmesh_assert(*returnval == *this);
#endif
    return returnval;
  }

  /**
   * Destructor.

  virtual dof_id_type n_nodes () const override final
  { return _n_nodes; }

  virtual dof_id_type parallel_n_nodes () const override final
  { return _n_nodes; }

  virtual dof_id_type max_node_id () const override final
  { return cast_int<dof_id_type>(_nodes.size()); }

  virtual dof_id_type n_elem () const override final
  { return _n_elem; }

  virtual dof_id_type parallel_n_elem () const override final
  { return _n_elem; }

  virtual dof_id_type n_active_elem () const override final;

  DECLARE_ELEM_ITERATORS(pid_, processor_id_type pid, pid)
  DECLARE_ELEM_ITERATORS(type_, ElemType type, type)

  DECLARE_ELEM_ITERATORS(active_subdomain_, subdomain_id_type sid, sid)
  DECLARE_ELEM_ITERATORS(active_subdomain_set_, std::set<subdomain_id_type> ss, ss)

  DECLARE_ELEM_ITERATORS(not_active_,,);

  DECLARE_ELEM_ITERATORS(active_pid_, processor_id_type pid, pid)
  DECLARE_ELEM_ITERATORS(local_level_, unsigned int level, level)
  DECLARE_ELEM_ITERATORS(local_not_level_, unsigned int level, level)
  DECLARE_ELEM_ITERATORS(active_local_subdomain_, subdomain_id_type sid, sid)
  DECLARE_ELEM_ITERATORS(active_local_subdomain_set_, std::set<subdomain_id_type> ss, ss)

  // Backwards compatibility
  virtual SimpleRange<element_iterator> active_subdomain_elements_ptr_range(subdomain_id_type sid) override final { return active_subdomain_element_ptr_range(sid); }
  virtual SimpleRange<const_element_iterator> active_subdomain_elements_ptr_range(subdomain_id_type sid) const override final { return active_subdomain_element_ptr_range(sid); }
  virtual SimpleRange<element_iterator> active_local_subdomain_elements_ptr_range(subdomain_id_type sid) override final { return active_local_subdomain_element_ptr_range(sid); }
  virtual SimpleRange<const_element_iterator> active_local_subdomain_elements_ptr_range(subdomain_id_type sid) const override final { return active_local_subdomain_element_ptr_range(sid); }
  virtual SimpleRange<element_iterator> active_subdomain_set_elements_ptr_range(std::set<subdomain_id_type> ss) override final { return active_subdomain_set_element_ptr_range(ss); }
  virtual SimpleRange<const_element_iterator> active_subdomain_set_elements_ptr_range(std::set<subdomain_id_type> ss) const override final { return active_subdomain_set_element_ptr_range(ss); }

void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T> & B)
{
  if (this->use_blas_lapack)
    this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
  else
    this->_right_multiply_transpose(B);
}

      LibmeshPetscCall(VecGhostRestoreLocalForm(_vec,&localrep));
    }
  else
    this->_type = SERIAL;
}

  if (this->comm().size() == 1)
    {
      libmesh_assert(ghost.empty());
      this->init(n, n_local, fast, ptype);
      return;
    }

  PetscInt petsc_n=static_cast<PetscInt>(n);

 * access, and the PETSc \p Mat returned by \p mat(), which is what solvers see.
 */
template <typename T>
class ReducedSparseMatrix : public PetscMatrixBase<T>
{
public:
  ReducedSparseMatrix(const Parallel::Communicator & comm_in, const ReducedDofMap & reduced_dof_map);

                    const numeric_index_type noz = 10,
                    const numeric_index_type blocksize = 1) override;

  virtual bool require_sparsity_pattern() const override { return false; }

  virtual void clear() noexcept override;

  /**
   * \returns The reduced degree of freedom map
   */
  const ReducedDofMap & reduced_dof_map() const { return _reduced_dof_map; }

  /**
   * \returns The wrapped matrix on the reduced space

 * that are used to map global system vector and matrix data to the condensed system vector and
 * matrix data
 */
class StaticCondensationDofMap : public ReducedDofMap
{
public:
  StaticCondensationDofMap(MeshBase & mesh, System & system, const DofMap & dof_map);

        continue;
      const Elem * parent = elem->parent();
      if (!parent || !elem->active())
        continue;

      ghost_objects_from_proc[obj_procid]++;
    }

        continue;
      const Elem * parent = elem->parent();
      if (!parent || !elem->active())
        continue;

      requested_objs_id[obj_procid].push_back(elem->id());
      requested_objs_parent_id_child_num[obj_procid].emplace_back

              query_id[i] = node.id();
            }
          else
            query_id[i] = DofObject::invalid_id;
        }

      // Gather whatever data the user wants

   * a negative value of the argument indicates we are
   * not performing a bounding box solve.
   */
  virtual void set_eigensolver_properties(int) {}

  /**
   * Set the name of the associated RB system --- we need

  /**
   * Get/set SCM_training_tolerance: tolerance for SCM greedy.
   */
  Real get_SCM_training_tolerance() const                         { return SCM_training_tolerance; }
  void set_SCM_training_tolerance(Real SCM_training_tolerance_in) { this->SCM_training_tolerance = SCM_training_tolerance_in; }

  /**
   * Perform the SCM greedy algorithm to develop a lower bound

   * default is does nothing. Override in subclass to attach a specific
   * vector.
   */
  virtual void attach_deflation_space() {}

protected:

   * indicator to be used in the SCM greedy.
   * Override in subclasses to specialize behavior.
   */
  virtual Real SCM_greedy_error_indicator(Real LB, Real UB) { return fabs(UB-LB)/fabs(UB); }

  //----------- PROTECTED DATA MEMBERS -----------//

   * Sets the position of the spectrum.
   */
  void set_position_of_spectrum (PositionOfSpectrum pos)
  {_position_of_spectrum= pos;}

  void set_position_of_spectrum (Real pos);
  void set_position_of_spectrum (Real pos, PositionOfSpectrum target);

  /**
   * \returns A constant reference to the system we are solving.
   */
  const sys_type & system () const { return _system; }

  /**
   * \returns A writable reference to the system we are solving.

  /**
   * Whether the solve eliminates constrained degrees of freedom
   */
  bool eliminating_constraints() const { return static_cast<IS>(_retained_is) != nullptr; }

  /**
   * Makes the system's \p current_local_solution the full iterate that \p x represents: \p x

#endif

private:
  friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
  friend void libmesh_petsc_snes_finish_residual (ResidualContext & rc,
                                                  Vec x,
                                                  Vec r,

  /**
   * \returns The number of converged eigenpairs.
   */
  unsigned int get_n_converged () const {return _n_converged_eigenpairs;}

  /**
   * \returns The number of eigen solver iterations.

                {
                  erase_nonhanging_vars(node, remaining_vars, vertices);
                  if (remaining_vars.empty())
                    continue;
                }

              erase_copied_vars(node, false, remaining_vars);

                {
                  erase_nonhanging_vars(node, remaining_vars, vertices);
                  if (remaining_vars.empty())
                    continue;
                }

              erase_copied_vars(node, false, remaining_vars);

   * @returns A prefix that may be applied to solver options. Note that this
   * prefix is only used if \p prefix_with_name()
   */
  std::string prefix() const { return this->name() + "_"; }

  /**
   * Request that static condensation be performed for this system

    Val & operator*() const { return it->second; }

    index_t index() const { return it->first; }

    veclike_iterator & operator++() { ++it; return *this; }

namespace libMesh
{
ConstraintEliminationDofMap::ConstraintEliminationDofMap(MeshBase & mesh,
                                                         System & system,
                                                         const DofMap & dof_map)
  : ReducedDofMap(mesh, system, dof_map, "constraint_eliminated")
{
  libmesh_experimental();
}

ConstraintEliminationDofMap::~ConstraintEliminationDofMap() = default;

bool ConstraintEliminationDofMap::retain_field_dof(const Elem &,
                                                   unsigned int,
                                                   unsigned int,
                                                   const dof_id_type field_dof) const
{
  return this->is_unconstrained(field_dof);
}

bool ConstraintEliminationDofMap::retain_scalar_dof(const dof_id_type scalar_dof) const
{
  return this->is_unconstrained(scalar_dof);
}

bool ConstraintEliminationDofMap::is_unconstrained(const dof_id_type dof) const
{
#ifdef LIBMESH_ENABLE_CONSTRAINTS
  return !_dof_map.is_constrained_dof(dof);
#else
  libmesh_ignore(dof);
  return true;

     implicit_neighbor_dofs,
     _need_full_sparsity_pattern,
     calculate_constrained,
     reduced);

  Threads::parallel_reduce (mesh.active_local_element_stored_range(), *sp);

  _matrices.clear();
  if (_sc)
    _sc->clear();
  if (_constraint_elimination)
    _constraint_elimination->clear();
}

    {
      // Static condensation reduces the space SNES solves on, which excludes constrained dofs when
      // they are eliminated
      _sc->set_parent(_constraint_elimination.get());
      _sc->reinit();
    }
}

void DofMap::create_constraint_elimination(MeshBase & mesh, System & sys)
{
  _constraint_elimination = std::make_unique<ConstraintEliminationDofMap>(mesh, sys, *this);
}

void DofMap::reinit_constraint_elimination()
{
  if (_constraint_elimination)
    _constraint_elimination->reinit();
}

unsigned int DofMap::add_variable(System & sys,

namespace libMesh
{
ReducedDofMap::ReducedDofMap(MeshBase & mesh,
                             System & system,
                             const DofMap & dof_map,
                             std::string reduced_system_name)
  : DofMapBase(dof_map.comm()),
    _mesh(mesh),
    _system(system),
    _dof_map(dof_map),
    _reduced_system_name(std::move(reduced_system_name)),
    _parent(nullptr),
    _initialized(false),
    _reduced_system(nullptr)
{
}

ReducedDofMap::~ReducedDofMap() = default;

void ReducedDofMap::retained_dof_added(const Elem &, dof_id_type) {}

void ReducedDofMap::unretained_dof_found(const Elem &, dof_id_type) {}

bool ReducedDofMap::retains_field_dof(const Elem & elem,
                                      const unsigned int node_num,
                                      const unsigned int var_num,
                                      const dof_id_type field_dof) const
{
  return (!_parent || _parent->retains_field_dof(elem, node_num, var_num, field_dof)) &&
         this->retain_field_dof(elem, node_num, var_num, field_dof);
}

bool ReducedDofMap::retains_scalar_dof(const dof_id_type scalar_dof) const
{
  return (!_parent || _parent->retains_scalar_dof(scalar_dof)) &&
         this->retain_scalar_dof(scalar_dof);
}

void ReducedDofMap::retain_dof(const dof_id_type full_dof)
{
  if (_dof_map.local_index(full_dof))
    _local_retained_set.insert(full_dof);
  else
    _nonlocal_retained[_dof_map.dof_owner(full_dof)].insert(full_dof);
}

void ReducedDofMap::retain_dof_and_constraining_dofs(const Elem * const elem,
                                                     const dof_id_type full_dof,
                                                     std::unordered_set<dof_id_type> & elem_seen,
                                                     std::vector<dof_id_type> & elem_var_dofs)
{
  if (!elem_seen.insert(full_dof).second)
    return;

  // A constraining dof belongs to this space only if it belongs to the parent space; one the
  // parent eliminates still passes its own constraining dofs on
  if (_parent && _parent->find_reduced_dof(full_dof) == DofObject::invalid_id)
    {
      this->retain_constraining_dofs(elem, full_dof, elem_seen, elem_var_dofs);
      return;
    }

  this->retain_dof(full_dof);
  elem_var_dofs.push_back(full_dof);
  this->retained_dof_added(*elem, full_dof);

  this->retain_constraining_dofs(elem, full_dof, elem_seen, elem_var_dofs);
}

void ReducedDofMap::retain_constraining_dofs(const Elem * const elem,
                                             const dof_id_type full_dof,
                                             std::unordered_set<dof_id_type> & elem_seen,
                                             std::vector<dof_id_type> & elem_var_dofs)
{
#ifdef LIBMESH_ENABLE_CONSTRAINTS
  const auto & constraints = _dof_map.get_dof_constraints();
  if (const auto it = constraints.find(full_dof); it != constraints.end())
    for (const auto & [constraining_dof, weight] : it->second)
      {
        libmesh_ignore(weight);
        // Constraining dofs may themselves be constrained
        this->retain_dof_and_constraining_dofs(elem, constraining_dof, elem_seen, elem_var_dofs);
      }
#else
  libmesh_ignore(elem, full_dof, elem_seen, elem_var_dofs);
#endif
}

void ReducedDofMap::reinit()
{
  if (this->initialized())
    this->clear();

  // SCALAR dofs are always retained. They are owned by the last rank
  for (const auto vg : make_range(_dof_map.n_variable_groups()))
    if (const auto & vg_description = _dof_map.variable_group(vg);
        vg_description.type().family == SCALAR)
      {
        std::vector<dof_id_type> scalar_dof_indices;
        for (const auto vg_vn : make_range(vg_description.n_variables()))
          {
            _dof_map.SCALAR_dof_indices(scalar_dof_indices, vg_description.number(vg_vn));
            for (const auto dof : scalar_dof_indices)
              if (this->retains_scalar_dof(dof))
                this->retain_dof(dof);
          }
      }

  std::vector<dof_id_type> elem_dofs; // only used to satisfy the dof_indices API
  std::vector<dof_id_type> * elem_var_dofs = nullptr;
  std::unordered_set<dof_id_type> elem_seen;
  std::set<unsigned int> full_vars_present;

  auto scalar_dofs_functor = [this, &elem_seen, &elem_var_dofs](
                                 const Elem & elem,
                                 std::vector<dof_id_type> & dof_indices,
                                 const std::vector<dof_id_type> & scalar_dof_indices) {
    dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
    for (const auto dof : scalar_dof_indices)
      if (this->retains_scalar_dof(dof))
        this->retain_dof_and_constraining_dofs(&elem, dof, elem_seen, *elem_var_dofs);
      else
        {
          if (!_dof_map.local_index(dof))
            _nonlocal_unretained.insert(dof);
          this->retain_constraining_dofs(&elem, dof, elem_seen, *elem_var_dofs);
        }
  };

  auto field_dofs_functor = [this, &elem_seen, &elem_var_dofs](const Elem & elem,
                                                              const unsigned int node_num,
                                                              const unsigned int var_num,
                                                              std::vector<dof_id_type> & dof_indices,
                                                              const dof_id_type field_dof) {
    dof_indices.push_back(field_dof);
    if (this->retains_field_dof(elem, node_num, var_num, field_dof))
      this->retain_dof_and_constraining_dofs(&elem, field_dof, elem_seen, *elem_var_dofs);
    else
      {
        // Only dofs of the parent space belong to this space's problem
        if (!_parent || _parent->retains_field_dof(elem, node_num, var_num, field_dof))
          this->unretained_dof_found(elem, field_dof);
        if (!_dof_map.local_index(field_dof))
          _nonlocal_unretained.insert(field_dof);
        this->retain_constraining_dofs(&elem, field_dof, elem_seen, *elem_var_dofs);
      }
  };

  for (const auto * const elem : _mesh.active_local_element_ptr_range())
    {
      auto & elem_dofs_by_var = _elem_reduced_dofs[elem->id()];
      elem_dofs_by_var.resize(_dof_map.n_variables());
      elem_seen.clear();

      const auto sub_id = elem->subdomain_id();
      for (const auto vg : make_range(_dof_map.n_variable_groups()))
        {
          const auto & var_group = _dof_map.variable_group(vg);
          if (!var_group.active_on_subdomain(sub_id))
            continue;
          for (const auto v : make_range(var_group.n_variables()))
            {
              const auto var_num = var_group.number(v);
              elem_var_dofs = &elem_dofs_by_var[var_num];
              _dof_map.dof_indices(
                  elem, elem_dofs, var_num, scalar_dofs_functor, field_dofs_functor, elem->p_level());
              if (!elem_var_dofs->empty())
                full_vars_present.insert(var_num);
            }
        }
    }

  //
  // A rank may own a retained dof that none of its local elements use, for instance a mid-edge
  // node of a coarse element that is a vertex of a finer element owned elsewhere. Every rank that
  // retains a nonlocal dof therefore tells the dof's owner
  //
  std::unordered_map<processor_id_type, std::vector<dof_id_type>> nonlocal_retained;
  for (const auto & [pid, set] : _nonlocal_retained)
    nonlocal_retained[pid].assign(set.begin(), set.end());
  _nonlocal_retained.clear();

  TIMPI::push_parallel_vector_data(
      this->comm(),
      nonlocal_retained,
      [this](processor_id_type, const std::vector<dof_id_type> & dofs) {
        _local_retained_set.insert(dofs.begin(), dofs.end());
      });

  _local_retained_dofs.assign(_local_retained_set.begin(), _local_retained_set.end());
  _local_retained_set.clear();
  std::sort(_local_retained_dofs.begin(), _local_retained_dofs.end());

  //
  // Number the reduced space
  //
  this->compute_dof_info(_local_retained_dofs.size());

  const auto full_first = _dof_map.first_dof();
  _local_full_to_reduced.assign(_dof_map.n_local_dofs(), DofObject::invalid_id);
  for (const auto i : index_range(_local_retained_dofs))
    _local_full_to_reduced[_local_retained_dofs[i] - full_first] = this->first_dof() + i;

  // Every nonlocal dof on a local element gets its reduced index, which is invalid_id for an
  // eliminated dof, so that find_reduced_dof() can tell eliminated dofs from unknown ones
  for (const auto dof : _nonlocal_unretained)
    nonlocal_retained[_dof_map.dof_owner(dof)].push_back(dof);
  _nonlocal_unretained.clear();

  auto gather_reduced = [this](processor_id_type,
                               const std::vector<dof_id_type> & full_dofs,
                               std::vector<dof_id_type> & reduced_dofs) {
    reduced_dofs.resize(full_dofs.size());
    for (const auto i : index_range(full_dofs))
      reduced_dofs[i] = this->find_reduced_dof(full_dofs[i]);
  };
  auto store_reduced = [this](processor_id_type,
                              const std::vector<dof_id_type> & full_dofs,
                              const std::vector<dof_id_type> & reduced_dofs) {
    for (const auto i : index_range(full_dofs))
      _nonlocal_full_to_reduced[full_dofs[i]] = reduced_dofs[i];
  };
  TIMPI::pull_parallel_vector_data(this->comm(),
                                   nonlocal_retained,
                                   gather_reduced,
                                   store_reduced,
                                   &DofObject::invalid_id);

  // The parent indices of the local retained dofs, and the ghosts of a vector on this space
  _local_retained_parent_dofs.resize(_local_retained_dofs.size());
  for (const auto i : index_range(_local_retained_dofs))
    _local_retained_parent_dofs[i] = this->parent_dof(_local_retained_dofs[i]);
  for (const auto & [full_dof, reduced] : _nonlocal_full_to_reduced)
    {
      libmesh_ignore(full_dof);
      if (reduced != DofObject::invalid_id)
        _send_list.push_back(reduced);
    }
  std::sort(_send_list.begin(), _send_list.end());

  //
  // Determine the variables with any dofs in the reduced space
  //
  this->comm().set_union(full_vars_present);
  _reduced_vars.reserve(full_vars_present.size());
  unsigned int first_local_number = 0;
  for (const auto full_var_num : full_vars_present)
    {
      const auto & full_var = _dof_map.variable(full_var_num);
      _reduced_vars.push_back(Variable{nullptr,
                                       full_var.name(),
                                       cast_int<unsigned int>(_reduced_vars.size()),
                                       first_local_number,
                                       full_var.type()});
      first_local_number += _reduced_vars.back().n_components(_mesh);
    }

  // Renumber the element dofs, keeping only the variables in the reduced space
  for (auto & [elem_id, elem_dofs_by_var] : _elem_reduced_dofs)
    {
      libmesh_ignore(elem_id);
      std::vector<std::vector<dof_id_type>> reduced_dofs_by_var;
      reduced_dofs_by_var.reserve(full_vars_present.size());
      for (const auto full_var_num : full_vars_present)
        {
          auto & var_dofs = reduced_dofs_by_var.emplace_back(
              std::move(elem_dofs_by_var[full_var_num]));
          for (auto & dof : var_dofs)
            dof = this->reduced_dof(dof);
        }
      elem_dofs_by_var = std::move(reduced_dofs_by_var);
    }

  // The sparsity pattern of the reduced space, with rows indexed by the local full dofs
  _reduced_sp = _dof_map.build_sparsity(_mesh, /*calculate_constrained=*/false, this);
  const auto & nnz = _reduced_sp->get_n_nz();
  const auto & noz = _reduced_sp->get_n_oz();
  libmesh_assert_equal_to(nnz.size(), noz.size());

  // Exact per-row preallocation counts. These are critical when SCALAR dofs introduce dense rows,
  // since the maximum row count would then allocate a dense matrix
  _reduced_nnz.resize(_local_retained_dofs.size());
  _reduced_noz.resize(_local_retained_dofs.size());
  for (const auto i : index_range(_local_retained_dofs))
    {
      const dof_id_type local_full_i = _local_retained_dofs[i] - full_first;
      libmesh_assert_less(local_full_i, nnz.size());
      _reduced_nnz[i] = nnz[local_full_i];
      _reduced_noz[i] = noz[local_full_i];
    }

  //
  // The stand-in system gives the reduced variables a system number. Nodes hold no reduced dofs,
  // since reduced dofs are queried by element only
  //
  const std::vector<unsigned int> nvpg(_reduced_vars.size(), 1);

  // add_system returns the existing system if there is one, so this is safe on every reinit
  _reduced_system = &_system.get_equation_systems().add_system("Basic", _reduced_system_name);
  if (!_reduced_system->is_initialized())
    {
      _reduced_system->init();
      for (auto * const nd : _mesh.active_node_ptr_range())
        {
          nd->set_n_vars_per_group(_reduced_system->number(), nvpg);
          for (const auto g : index_range(nvpg))
            nd->set_n_comp_group(_reduced_system->number(), g, 0);
        }
    }

  // The stand-in system holds no solution data worth writing
  _reduced_system->hide_output() = true;

  _initialized = true;
}

dof_id_type ReducedDofMap::find_reduced_dof(const dof_id_type full_dof) const
{
  if (_dof_map.local_index(full_dof))
    return _local_full_to_reduced[full_dof - _dof_map.first_dof()];
  return libmesh_map_find(_nonlocal_full_to_reduced, full_dof);
}

dof_id_type ReducedDofMap::parent_dof(const dof_id_type full_dof) const
{
  return _parent ? _parent->find_reduced_dof(full_dof) : full_dof;
}

dof_id_type ReducedDofMap::parent_n_dofs() const
{
  return _parent ? _parent->n_dofs() : _dof_map.n_dofs();
}

dof_id_type ReducedDofMap::parent_n_local_dofs() const
{
  return _parent ? _parent->n_local_dofs() : _dof_map.n_local_dofs();
}

dof_id_type ReducedDofMap::parent_first_dof() const
{
  return _parent ? _parent->first_dof() : _dof_map.first_dof();
}

dof_id_type ReducedDofMap::parent_end_dof() const
{
  return _parent ? _parent->end_dof() : _dof_map.end_dof();
}

const std::vector<dof_id_type> & ReducedDofMap::parent_send_list() const
{
  return _parent ? _parent->send_list() : _dof_map.get_send_list();
}

dof_id_type ReducedDofMap::reduced_dof(const dof_id_type full_dof) const
{
  const auto reduced = this->find_reduced_dof(full_dof);
  libmesh_assert_msg(reduced != DofObject::invalid_id,
                     "Full dof " << full_dof << " is not retained in the reduced space");
  return reduced;
}

const std::vector<std::vector<dof_id_type>> &
ReducedDofMap::elem_dof_indices(const dof_id_type elem_id) const
{
  return libmesh_map_find(_elem_reduced_dofs, elem_id);
}

unsigned int ReducedDofMap::n_variables() const { return _reduced_vars.size(); }

const Variable & ReducedDofMap::variable(const unsigned int c) const
{
  libmesh_assert_less(c, _reduced_vars.size());
  return _reduced_vars[c];
}

void ReducedDofMap::dof_indices(const Elem * const elem,
                                std::vector<dof_id_type> & di,
                                const unsigned int vn,
                                int /*p_level*/) const
{
  di = libmesh_map_find(_elem_reduced_dofs, elem->id())[vn];
}

void ReducedDofMap::dof_indices(const Node * const, std::vector<dof_id_type> &, const unsigned int) const
{
  libmesh_error_msg("ReducedDofMap dof indices are only meant to be queried with elements, not "
                    "nodes");
}

void ReducedDofMap::clear()
{
  DofMapBase::clear();
  _local_retained_dofs.clear();
  _local_retained_parent_dofs.clear();
  _send_list.clear();
  _local_full_to_reduced.clear();
  _nonlocal_full_to_reduced.clear();
  _local_retained_set.clear();
  _nonlocal_retained.clear();
  _nonlocal_unretained.clear();
  _elem_reduced_dofs.clear();
  _reduced_vars.clear();
  _reduced_sp = nullptr;
  _reduced_nnz.clear();
  _reduced_noz.clear();
  _initialized = false;
}
}

              const bool implicit_neighbor_dofs_in,
              const bool need_full_sparsity_pattern_in,
              const bool calculate_constrained_in,
              const ReducedDofMap * const reduced_in) :
  ParallelObject(dof_map_in),
  dof_map(dof_map_in),
  dof_coupling(dof_coupling_in),

  implicit_neighbor_dofs(implicit_neighbor_dofs_in),
  need_full_sparsity_pattern(need_full_sparsity_pattern_in),
  calculate_constrained(calculate_constrained_in),
  reduced(reduced_in),
  sparsity_pattern(),
  nonlocal_pattern(),
  n_nz(),

  implicit_neighbor_dofs(other.implicit_neighbor_dofs),
  need_full_sparsity_pattern(other.need_full_sparsity_pattern),
  calculate_constrained(other.calculate_constrained),
  reduced(other.reduced),
  hashed_dof_sets(other.hashed_dof_sets),
  sparsity_pattern(),
  nonlocal_pattern(),

                                  std::vector<dof_id_type> & dofs_vi,
                                  unsigned int vi)
{
  if (this->reduced)
  {
    // We build a sparsity pattern that will match the size of the reduced system. This is so that
    // we have the data necessary to init the reduced system matrix

    // Unretained dofs still couple their constraining dofs, which constrained element matrices
    // and vectors act on
    dummy_unretained.clear();

    auto total_and_retained_from_scalar_dofs_functor =
        [&dofs_vi, this](const Elem & /*elem*/,
                         std::vector<dof_id_type> & dof_indices,
                         const std::vector<dof_id_type> & scalar_dof_indices)
    {
      dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
      for (const auto scalar_dof : scalar_dof_indices)
        (this->reduced->retains_scalar_dof(scalar_dof) ? dofs_vi : dummy_unretained)
            .push_back(scalar_dof);
    };

    auto total_and_retained_from_field_dofs_functor =

                        const dof_id_type field_dof)
    {
      dof_indices.push_back(field_dof);
      (this->reduced->retains_field_dof(functor_elem, node_num, var_num, field_dof)
           ? dofs_vi
           : dummy_unretained)
          .push_back(field_dof);
    };

    dof_map.dof_indices(elem,

#ifdef LIBMESH_ENABLE_CONSTRAINTS
    // Constraint rows reference unconstrained dofs only, so one level suffices
    const auto & constraints = dof_map.get_dof_constraints();
    for (const auto unretained : dummy_unretained)
      if (const auto it = constraints.find(unretained); it != constraints.end())
        for (const auto & [constraining_dof, weight] : it->second)
          {
            libmesh_ignore(weight);
            dofs_vi.push_back(constraining_dof);
          }
#endif
  }

      // want to handle neighbor links from subactive descendents.
      const unsigned int member_level = neigh_family_member->level();
      if (member_level < this_level)
        continue;

      // Ideally, the neighbor link ought to either be correct
      // already or ought to be to remote_elem.

          if (this->ancestor())
            continue;

          if (member_subactive && this->has_children())
            continue;
        }

      // If neigh is at a coarser level than us, some of neigh's

        {
          libmesh_assert(member_subactive);

          this->side_ptr(my_side, n);
          neigh_family_member->side_ptr(neigh_side, nn);

          if (*my_side != *neigh_side)
            continue;
        }

      neigh_family_member->set_neighbor(nn, this);

              // neighbor
              Elem * my_ancestor = this->parent();
              libmesh_assert(my_ancestor);
              while (!neigh->has_neighbor(my_ancestor))
                {
                  my_ancestor = my_ancestor->parent();
                  libmesh_assert(my_ancestor);

              // My neighbor may have descendants which consider me a
              // neighbor
              std::vector<Elem *> family;
              neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);

              for (auto & n : family)
                {
                  libmesh_assert (n);
                  if (n->is_remote())
                    continue;
                  unsigned int my_s = n->which_neighbor_am_i(this);
                  libmesh_assert_less (my_s, n->n_neighbors());
                  libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
                  // TODO: we may want to make remote_elem non-const.
                  n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
                }
            }
#endif

void Elem::total_family_tree_by_subneighbor (std::vector<Elem *> & family,
                                             Elem * neighbor,
                                             Elem * subneighbor,
                                             bool reset)
{
  ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
}

  std::unique_ptr<PointLocatorBase> point_locator;
  if (check_periodic_bcs || check_disjoint_bcs)
      point_locator = _mesh->sub_point_locator();

  std::set<const Elem *> periodic_elems_examined;
  const BoundaryInfo & binfo = _mesh->get_boundary_info();

      if (check_disjoint_bcs)
        {
          // Also ghost their disjoint neighbors
          for (auto s : elem->side_index_range())
          {
            for (const auto & [id, boundary_ptr] : *db)
            {
              if (!_mesh->get_boundary_info().has_boundary_id(elem, s, id))
                continue;

              unsigned int neigh_side = invalid_uint;
              const Elem * neigh =
                  db->neighbor(id, *point_locator, elem, s, &neigh_side);

              if (!neigh || neigh == remote_elem)
                continue;

              if (neigh->processor_id() != p)
                coupled_elements.emplace(neigh, nullcm);
            }
          }
        }

                    other_boundary_info._ns_id_to_name) ||
      !compare_maps(_es_id_to_name,
                    other_boundary_info._es_id_to_name))
    return false;

  return true;
}

  std::set<boundary_id_type> request_boundary_ids(_boundary_ids);
  request_boundary_ids.insert(invalid_id);
  if (!_mesh->is_serial())
    this->comm().set_union(request_boundary_ids);

  this->sync(request_boundary_ids,
             boundary_mesh);

   * easily misused.
   */
  if (!_mesh->is_serial())
    boundary_mesh.delete_remote_elements();

  /**
   * If the boundary_mesh is still serial, that means we *can't*

   */
  std::unique_ptr<MeshBase> mesh_copy;
  if (boundary_mesh.is_serial() && !_mesh->is_serial())
    mesh_copy = _mesh->clone();

  auto serializer = std::make_unique<MeshSerializer>
    (const_cast<MeshBase &>(*_mesh), boundary_mesh.is_serial());

  // leave those pointers dangling.  Fix them up if needed.
  if (mesh_copy.get())
    {
      for (auto & new_elem : boundary_mesh.element_ptr_range())
        {
          const dof_id_type interior_parent_id =
            new_elem->interior_parent()->id();

          if (!mesh_copy->query_elem_ptr(interior_parent_id))
            new_elem->set_interior_parent
              (const_cast<RemoteElem *>(remote_elem));
        }
    }

  // Don't repartition this mesh; we want it to stay in sync with the

      if (elem->has_children())
        for (auto c : make_range(elem->n_children()))
          if (elem->child_ptr(c) == remote_elem &&
              elem->is_child_on_side(c, s))
            {
              for (auto sc : make_range(new_elem->n_children()))
                if (!new_elem->raw_child_ptr(sc))
                  new_elem->add_child
                    (const_cast<RemoteElem*>(remote_elem), sc);
            }
#endif

      // On non-local elements on DistributedMesh we might have
      // RemoteElem neighbor links to construct
      if (!_mesh->is_serial() &&
          (elem->processor_id() != this->processor_id()))
        {
          const unsigned short n_nodes = elem->n_nodes();

          const unsigned short bdy_n_sides = new_elem->n_sides();
          const unsigned short bdy_n_nodes = new_elem->n_nodes();

          // Check every interior side for a RemoteElem
          for (auto interior_side : elem->side_index_range())
            {
              // Might this interior side have a RemoteElem that
              // needs a corresponding Remote on a boundary side?
              if (elem->neighbor_ptr(interior_side) != remote_elem)
                continue;

              // Which boundary side?
              for (unsigned short boundary_side = 0;
                   boundary_side != bdy_n_sides; ++boundary_side)
                {
                  // Look for matching node points.  This is safe in
                  // *this* context.
                  bool found_all_nodes = true;
                  for (unsigned short boundary_node = 0;
                       boundary_node != bdy_n_nodes; ++boundary_node)
                    {
                      if (!new_elem->is_node_on_side(boundary_node,
                                                     boundary_side))
                        continue;

                      bool found_this_node = false;
                      for (unsigned short interior_node = 0;
                           interior_node != n_nodes; ++interior_node)
                        {
                          if (!elem->is_node_on_side(interior_node,
                                                     interior_side))
                            continue;

                          if (new_elem->point(boundary_node) ==
                              elem->point(interior_node))

                              break;
                            }
                        }
                      if (!found_this_node)
                        {
                          found_all_nodes = false;
                          break;
                        }
                    }

                  if (found_all_nodes)
                    {
                      new_elem->set_neighbor
                        (boundary_side,
                         const_cast<RemoteElem *>(remote_elem));
                      break;
                    }

  // Don't add the same ID twice
  for (const auto & pr : as_range(_boundary_shellface_id.equal_range(elem)))
    if (pr.second.first == shellface &&
        pr.second.second == id)
      return;

  _boundary_shellface_id.emplace(elem, std::make_pair(shellface, id));

#ifdef LIBMESH_ENABLE_AMR

          while (p != nullptr)
          {
            const Elem * parent = p->parent();
            if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side))
              break;
            p = parent;
          }
#endif
          // We're on that side of our top_parent; return it
          if (!p)
            returnval.push_back(side);
        }
        // Otherwise we trust what we got and return the side
        else

  std::size_t n_edge_bcs=0;

  for (const auto & pr : _boundary_edge_id)
    if (pr.first->processor_id() == this->processor_id())
      n_edge_bcs++;

  this->comm().sum (n_edge_bcs);

  std::size_t n_shellface_bcs=0;

  for (const auto & pr : _boundary_shellface_id)
    if (pr.first->processor_id() == this->processor_id())
      n_shellface_bcs++;

  this->comm().sum (n_shellface_bcs);

              if (!mesh_is_serial)
                {
                  const processor_id_type proc_id =
                    side->node_ptr(i)->processor_id();
                  if (proc_id != my_proc_id)
                    nodes_to_push[proc_id].emplace(side->node_id(i), bcid);
                }
            }
        }

  // Otherwise we need to push ghost node bcids to their owners, then
  // pull ghost node bcids from their owners.

  for (auto & [proc_id, s] : nodes_to_push)
    {
      node_vecs_to_push[proc_id].assign(s.begin(), s.end());
      s.clear();
    }

  auto nodes_action_functor =
    [this]
    (processor_id_type,
     const vec_type & received_nodes)
    {
      for (const auto & [dof_id, bndry_id] : received_nodes)
        this->add_node(_mesh->node_ptr(dof_id), bndry_id);
    };

  Parallel::push_parallel_vector_data
    (this->comm(), node_vecs_to_push, nodes_action_functor);

  // At this point we should know all the BCs for our own nodes; now
  // we need BCs for ghost nodes.

    node_ids_requested;

  // Determine what nodes we need to request
  for (const auto & node : _mesh->node_ptr_range())
    {
      const processor_id_type pid = node->processor_id();
      if (pid != my_proc_id)
        node_ids_requested[pid].push_back(node->id());
    }

  typedef std::vector<boundary_id_type> datum_type;

  auto node_bcid_gather_functor =
    [this]
    (processor_id_type,
     const std::vector<dof_id_type> & ids,
     std::vector<datum_type> & data)
    {
      const std::size_t query_size = ids.size();
      data.resize(query_size);

      for (std::size_t i=0; i != query_size; ++i)
        this->boundary_ids(_mesh->node_ptr(ids[i]), data[i]);
    };

  auto node_bcid_action_functor =
    [this]
    (processor_id_type,
     const std::vector<dof_id_type> & ids,
     const std::vector<datum_type> & data)
    {
      for (auto i : index_range(ids))
        this->add_node(_mesh->node_ptr(ids[i]), data[i]);
    };

  datum_type * datum_type_ex = nullptr;
  Parallel::pull_parallel_vector_data
    (this->comm(), node_ids_requested, node_bcid_gather_functor,
     node_bcid_action_functor, datum_type_ex);
}

        // subdomain id.
        if (neigh == remote_elem)
          {
            const processor_id_type pid = elem->processor_id();
            sides_requested[pid].emplace_back(elem->id(), side);
            requested_boundary_ids[pid].push_back(std::move(matched_ids));
            continue;
          }

        add_matched_side(elem, side, matched_ids);
      } // end for side

  // below rather than skipping it when its own request list is empty.

  auto neighbor_subdomain_gather_functor =
    [this]
    (processor_id_type,
     const std::vector<elem_side_type> & elem_sides,
     std::vector<subdomain_id_type> & neighbor_subdomains)
    {
      neighbor_subdomains.resize(elem_sides.size());
      for (auto i : index_range(elem_sides))
        {
          const Elem & elem = _mesh->elem_ref(elem_sides[i].first);
          const Elem * neigh = elem.neighbor_ptr(elem_sides[i].second);

          // The requester saw a RemoteElem neighbor here, so there is
          // a neighbor, and we own elem, so our connectivity for it is

          libmesh_assert(neigh);
          libmesh_assert(neigh != remote_elem);

          neighbor_subdomains[i] = neigh->subdomain_id();
        }
    };

  auto neighbor_subdomain_action_functor =
    [this, &requested_boundary_ids, &add_matched_side]
    (processor_id_type pid,
     const std::vector<elem_side_type> & elem_sides,
     const std::vector<subdomain_id_type> & neighbor_subdomains)
    {
      const auto & bids = libmesh_map_find(requested_boundary_ids, pid);
      libmesh_assert_equal_to(bids.size(), elem_sides.size());

      for (auto i : index_range(elem_sides))
        {
          const Elem * elem = _mesh->elem_ptr(elem_sides[i].first);

          // If the side turned out to be truly interior after all,
          // skip it.
          if (neighbor_subdomains[i] != elem->subdomain_id())
            add_matched_side(elem, elem_sides[i].second, bids[i]);
        }
    };

}


void DistributedMesh::own_node (Node & n)
{
  // This had better be a node in our mesh
  libmesh_assert(_nodes[n.id()] == &n);

  _nodes[n.id()] = nullptr;
  _n_nodes--;

  n.set_id(DofObject::invalid_id);
  n.processor_id() = this->processor_id();

  this->add_node(&n);
}


Node * DistributedMesh::add_node (Node * n)

              // else does, better figure out who should own it next.
              if (!used_nodes.count(n))
                {
                  if (auto it = repartitioned_node_pids.find(n);
                      sender_could_become_owner)
                    {
                      if (it != repartitioned_node_pids.end() &&
                          pid < it->second)
                        it->second = pid;
                      else
                        repartitioned_node_pids[n] = pid;
                    }
                  else
                    if (it == repartitioned_node_pids.end())
                      repartitioned_node_pids[n] =
                        DofObject::invalid_processor_id;

                  repartitioned_node_sets_to_push[pid].insert(n);
                }
            }
        };

      // Repartition (what used to be) our own nodes first
      for (auto & [n, p] : repartitioned_node_pids)
        {
          Node & node = this->node_ref(n);
          libmesh_assert_equal_to(node.processor_id(), this->processor_id());
          libmesh_assert_not_equal_to_msg(p, DofObject::invalid_processor_id, "Node " << n << " is lost?");
          node.processor_id() = p;
        }

      // Then push to repartition others' ghosted copies.

      for (auto & [p, nodeset] : repartitioned_node_sets_to_push)
        {
          auto & rn_vec = repartitioned_node_vecs[p];
          for (auto n : nodeset)
            rn_vec.emplace_back(n, repartitioned_node_pids[n]);
        }

      auto repartition_node_functor =
        [this]
        (processor_id_type libmesh_dbg_var(pid),
         const std::vector<std::pair<dof_id_type, processor_id_type>> & ids_and_pids)
        {
          for (auto [n, p] : ids_and_pids)
            {
              libmesh_assert_not_equal_to(p, DofObject::invalid_processor_id);
              Node & node = this->node_ref(n);
              libmesh_assert_equal_to(node.processor_id(), pid);
              node.processor_id() = p;
            }
        };

void DistributedMesh::fix_broken_node_and_element_numbering ()
{
  // We can't use range-for here because we need access to the special
  // iterators' methods, not just to their dereferenced values.

  // Nodes first
  for (auto pr = this->_nodes.begin(),
           end = this->_nodes.end(); pr != end; ++pr)
    {
      Node * n = *pr;
      if (n != nullptr)
        {
          const dof_id_type id = pr.index();
          n->set_id() = id;
          libmesh_assert_equal_to(this->node_ptr(id), n);
        }
    }

  // Elements next
  for (auto pr = this->_elements.begin(),
           end = this->_elements.end(); pr != end; ++pr)
    {
      Elem * e = *pr;
      if (e != nullptr)
        {
          const dof_id_type id = pr.index();
          e->set_id() = id;
          libmesh_assert_equal_to(this->elem_ptr(id), e);
        }
    }
}

  // Do nothing if no variables are detected
  if (count == 0)
    return;

  // Second read the actual names and convert them into a format we can use
  NamesData names_table(count, libmesh_max_str_length);

      const auto & constraint_rows = mesh->get_constraint_rows();

      std::unordered_set<const Elem *> constraining_nodes_elems;
      for (const Elem * elem : connected_elements)
        {
          for (const Node & node : elem->node_ref_range())
            {
              // Retain all elements containing constraining nodes
              if (const auto it = constraint_rows.find(&node);
                  it != constraint_rows.end())
                for (auto & p : it->second)
                  {
                    const Elem * constraining_elem = p.first.first;
                    libmesh_assert(constraining_elem ==
                                   mesh->elem_ptr(constraining_elem->id()));
                    if (!connected_elements.count(constraining_elem) &&

            }
        }

      newer_connected_elements.insert(constraining_nodes_elems.begin(),
                                      constraining_nodes_elems.end());
    }

          if (e->active() && e->has_children())
            {
              std::vector<const Elem *> subactive_family;
              e->total_family_tree(subactive_family);
              for (const auto & f : subactive_family)
                {
                  libmesh_assert(f != remote_elem);
                  if (!connected_elements.count(f) &&
                      !new_connected_elements.count(f))
                    newer_connected_elements.insert(f);
                }
            }
        };

      for (auto & [node, row] : constraint_rows)
        {
          if (!connected_nodes.count(node))
            continue;

          serialized_row_type serialized_row;
          for (auto [elem_and_node, coef] : row)
            serialized_row.emplace_back(std::make_pair(elem_and_node.first->id(),
                                                       elem_and_node.second),
                                        coef);

          all_constraint_rows_to_send[pid].emplace_back
            (node->id(), std::move(serialized_row));
        }
    }

  if (have_constraint_rows)
    {
      auto constraint_row_action =
        [&mesh, &constraint_rows]
        (processor_id_type /* src_pid */,
         const std::vector<std::pair<dof_id_type, serialized_row_type>> rows)
        {
          for (auto & [node_id, serialized_row] : rows)
            {
              MeshBase::constraint_rows_mapped_type row;
              for (auto [elem_and_node, coef] : serialized_row)
                row.emplace_back(std::make_pair(mesh.elem_ptr(elem_and_node.first),
                                                elem_and_node.second),
                                 coef);

              constraint_rows[mesh.node_ptr(node_id)] = row;
            }
        };

      TIMPI::push_parallel_vector_data(mesh.comm(),
                                       all_constraint_rows_to_send,
                                       constraint_row_action);

      // the parent's owner needs us to send them.
      const processor_id_type their_proc_id = elem->parent()->processor_id();
      if (their_proc_id != proc_id)
        coarsening_elements_to_ghost[their_proc_id].push_back(elem);
    }

  std::map<processor_id_type, std::vector<const Node *>> all_nodes_to_send;

          // Make some fake element iterators defining this vector of
          // elements
          Elem * const * elempp = const_cast<Elem * const *>(elems.data());
          Elem * const * elemend = elempp+elems.size();
          const MeshBase::const_element_iterator elem_it =
            MeshBase::const_element_iterator(elempp, elemend, Predicates::NotNull<Elem * const *>());
          const MeshBase::const_element_iterator elem_end =
            MeshBase::const_element_iterator(elemend, elemend, Predicates::NotNull<Elem * const *>());

          for (auto & gf : as_range(mesh.ghosting_functors_begin(),
                                    mesh.ghosting_functors_end()))
            {
              GhostingFunctor::map_type elements_to_ghost;
              libmesh_assert(gf);
              (*gf)(elem_it, elem_end, p, elements_to_ghost);

              // We can ignore the CouplingMatrix in ->second, but we
              // need to ghost all the elements in ->first.
              for (auto & pr : elements_to_ghost)
                {
                  const Elem * elem = pr.first;
                  libmesh_assert(elem);
                  while (elem)
                    {
                      libmesh_assert(elem != remote_elem);
                      elements_to_send.insert(elem);
                      for (auto & n : elem->node_ref_range())
                        nodes_to_send.insert(&n);
                      elem = elem->parent();
                    }
                }
            }

          all_nodes_to_send[p].assign(nodes_to_send.begin(), nodes_to_send.end());
          all_elems_to_send[p].assign(elements_to_send.begin(), elements_to_send.end());
        }
    }

               std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
        serialized_rows;

      for (auto & row : constraint_rows)
        {
          const Node * node = row.first;
          const dof_id_type rowid = node->id();
          libmesh_assert(node == mesh.node_ptr(rowid));

          std::vector<std::tuple<dof_id_type, unsigned int, Real>>
            serialized_row;
          for (auto & entry : row.second)
            serialized_row.push_back
              (std::make_tuple(entry.first.first->id(),
                               entry.first.second, entry.second));

          serialized_rows.emplace(rowid, std::move(serialized_row));
        }

      if (root_id == DofObject::invalid_processor_id)
        mesh.comm().set_union(serialized_rows);
      else
        mesh.comm().set_union(serialized_rows, root_id);

      if (root_id == DofObject::invalid_processor_id ||
          root_id == mesh.processor_id())
        {
          for (auto & row : serialized_rows)
            {
              const dof_id_type rowid = row.first;
              const Node * node = mesh.node_ptr(rowid);

              std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
                deserialized_row;
              for (auto & entry : row.second)
                deserialized_row.push_back
                  (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
                                                 std::get<1>(entry)),
                                                 std::get<2>(entry)));

              constraint_rows.emplace(node, deserialized_row);
            }
        }
#ifdef DEBUG

          {
            // But we might have gotten a definitive id from a
            // different request
            if (!definitive_ids.count(mesh.node_ptr(old_id)))
              data_changed = true;
          }
        else

{
  typedef std::pair<unsigned char,unsigned char> datum;

  SyncPLevels(MeshBase & _mesh) :
    mesh(_mesh) {}

  MeshBase & mesh;

  // Find the p_level of each requested Elem
  void gather_data (const std::vector<dof_id_type> & ids,
                    std::vector<datum> & ids_out) const
  {
    ids_out.reserve(ids.size());

    for (const auto & id : ids)
      {
        Elem & elem = mesh.elem_ref(id);
        ids_out.push_back
          (std::make_pair(cast_int<unsigned char>(elem.p_level()),
                          static_cast<unsigned char>(elem.p_refinement_flag())));
      }
  }

  void act_on_data (const std::vector<dof_id_type> & old_ids,
                    const std::vector<datum> & new_p_levels) const
  {
    for (auto i : index_range(old_ids))
      {
        Elem & elem = mesh.elem_ref(old_ids[i]);
        // Make sure these are consistent
        elem.hack_p_level_and_refinement_flag
          (new_p_levels[i].first,
           static_cast<Elem::RefinementState>(new_p_levels[i].second));
        // Make sure parents' levels are consistent
        elem.set_p_level(new_p_levels[i].first);
      }
  }
};
#endif // LIBMESH_ENABLE_AMR

// ------------------------------------------------------------
#ifdef LIBMESH_ENABLE_AMR
void MeshCommunication::make_p_levels_parallel_consistent(MeshBase & mesh)
{
  // This function must be run on all processors at once
  libmesh_parallel_only(mesh.comm());

  SyncPLevels syncplevels(mesh);
  Parallel::sync_dofobject_data_by_id
    (mesh.comm(), mesh.elements_begin(), mesh.elements_end(),
     syncplevels);
}
#endif // LIBMESH_ENABLE_AMR

        if (node.processor_id() > proc_ids[i])
          {
            data_changed = true;
            node.processor_id() = proc_ids[i];
          }
      }

  // If cross_section is distributed, so is its extrusion
  if (!cross_section.is_serial())
    mesh.delete_remote_elements();

  // We know a priori how many elements we'll need
  mesh.reserve_elem(nz*orig_elem);

      Elem * elem = mesh.query_elem_ptr(e_id);

      if (!elem)
        continue;

      const unsigned int max_permutation = elem->n_permutations();
      if (!max_permutation)

              // final tet.
              if (next_subelem == 3)
                {
                  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
                  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
                  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
                  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
                  subelem[next_subelem]->orient(&boundary_info);
                  ++next_subelem;

                  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
                  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
                  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
                  subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
                  subelem[next_subelem]->orient(&boundary_info);
                  ++next_subelem;

                  // We don't need the 6th tet after all

      std::array<double, 3> point_val;

      // We should only be getting new nodes if we asked for them
      if (!_desired_volume)
        {
          std::cout <<
            "NetGen output " << n_points <<

          libmesh_error();
        }
      else
        for (auto i : make_range(old_nodes, n_points))
          {
            // i+1 since ng uses ONE-BASED numbering
            Ng_GetPoint (ngmesh, i+1, point_val.data());
            const Point p(point_val[0], point_val[1], point_val[2]);
            Node * n_new = this->_mesh.add_point(p);
            const dof_id_type n_new_id = n_new->id();
            ng_to_libmesh_id[i+1] = n_new_id;
          }
    }

              // our change has to propagate to neighboring
              // processors.
              if (my_flag_changed && !_mesh.is_serial())
                for (auto n : elem->side_index_range())
                  {
                    Elem * neigh =
                      topological_neighbor(elem, point_locator.get(), n);

                    if (!neigh)
                      continue;
                    if (neigh == remote_elem ||
                        neigh->processor_id() !=
                        this->processor_id())
                      {
                        compatible_with_refinement = false;
                        break;
                      }
                    // FIXME - for non-level one meshes we should
                    // test all descendants
                    if (neigh->has_children())
                      for (auto & child : neigh->child_ref_range())
                        if (&child == remote_elem ||
                            child.processor_id() !=
                            this->processor_id())
                          {
                            compatible_with_refinement = false;
                            break;
                          }
                  }
            }

                break;
              if (child.processor_id() != elem->processor_id())
                {
                  uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
                  break;
                }
            }

          for (const auto & id : my_uncoarsenable_parents)
            {
              Elem & elem = _mesh.elem_ref(id);
              libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
                             elem.refinement_flag() == Elem::COARSEN_INACTIVE);
              elem.set_refinement_flag(Elem::INACTIVE);

  // levels changed in sync
  if (mesh_p_changed && !_mesh.is_serial())
    {
      MeshCommunication().make_p_levels_parallel_consistent (_mesh);
    }

  return (mesh_changed || mesh_p_changed);

          if (elem->p_refinement_flag() == Elem::REFINE &&
              elem->active())
            {
              elem->set_p_level(elem->p_level()+1);
              elem->set_p_refinement_flag(Elem::JUST_REFINED);
              mesh_p_changed = true;
            }
        }
    }

  if (mesh_p_changed && !_mesh.is_replicated())
    {
      MeshCommunication().make_p_levels_parallel_consistent (_mesh);
    }

  // Clear the _new_nodes_map and _unused_elements data structures.

            // independent of element traversal.
            if (!mesh.is_replicated())
              {
                side_elem->set_id(max_orig_id + max_sides*elem->id() + s);
#ifdef LIBMESH_ENABLE_UNIQUE_ID
                side_elem->set_unique_id(max_unique_id + max_sides*elem->id() + s);
#endif
              }
          }

        typedef std::unordered_multimap<key_type, val_type> map_type;
        map_type side_to_elem_map;

        std::unique_ptr<Elem> my_side, their_side;

        for (auto & elem : elems_to_add)
          {
            for (auto s : elem->side_index_range())
              {
                if (elem->neighbor_ptr(s))
                  continue;
                const dof_id_type key = elem->low_order_key(s);
                auto bounds = side_to_elem_map.equal_range(key);
                if (bounds.first != bounds.second)
                  {
                    elem->side_ptr(my_side, s);
                    while (bounds.first != bounds.second)

        // At this point we *should* have a match for everything, so
        // anything we don't have a match for is remote.
        for (auto & elem : elems_to_add)
          for (auto s : elem->side_index_range())
            if (!elem->neighbor_ptr(s))
              elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
      }

    // Remove volume and edge elements
    for (Elem * elem : elems_to_delete)

    {
      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

  processor_id_type pid_broadcasting_names = this->processor_id();
  const std::size_t n_names = result.size();
  if (!n_names)
    pid_broadcasting_names = DofObject::invalid_processor_id;

  libmesh_assert(this->comm().semiverify
                 (n_names ? nullptr : &n_names));

      if (const processor_id_type pid = elem.processor_id();
          pid != _mesh.processor_id() && !_mesh.is_serial() &&
          !_desired_volume_func.get() &&
          elem.id() == coarse_id)
        edge_queries[pid].emplace_back(vertices.first->id(),
                                       vertices.second->id());

      // Test should_refine, but also test for any edges that were
      // already refined by neighbors, which might happen even if

  // some unsplit edges that could have been split by their
  // remote_elem neighbors, and we'll need to query those.
  auto edge_gather_functor =
    [this]
    (processor_id_type,
     const std::vector<std::pair<dof_id_type, dof_id_type>> & edges,
     std::vector<refinement_datum> & edge_refinements)
    {
      // Fill those requests
      const std::size_t query_size = edges.size();
      edge_refinements.resize(query_size);
      for (std::size_t i=0; i != query_size; ++i)
        {
          auto vertex_ids = edges[i];
          std::pair<Node *, Node *> vertices
            {_mesh.node_ptr(vertex_ids.first),
             _mesh.node_ptr(vertex_ids.second)};
          fill_refinement_datum(vertices, edge_refinements[i]);
        }
    };

  auto edge_action_functor =
    [this]
    (processor_id_type,
     const std::vector<std::pair<dof_id_type, dof_id_type>> &,
     const std::vector<refinement_datum> & edge_refinements
    )
    {
      for (const auto & one_edges_refinements : edge_refinements)
        for (const auto & refinement : one_edges_refinements)
          {
            std::pair<Node *, Node *> vertices
              {_mesh.node_ptr(std::get<0>(refinement)),
               _mesh.node_ptr(std::get<1>(refinement))};

            if ((Point&)(*vertices.first) >
                (Point&)(*vertices.second))
              std::swap(vertices.first, vertices.second);

            if (auto it = new_nodes.find(vertices);
                it != new_nodes.end())
              {
                _mesh.renumber_node(it->second->id(),
                                    std::get<2>(refinement));
                it->second->processor_id() = std::get<3>(refinement);
              }
            else
              {
                Node * new_node =
                  _mesh.add_point((*(Point*)vertices.first +
                                   *(Point*)vertices.second)/2,
                                  std::get<2>(refinement),
                                  std::get<3>(refinement));
                new_nodes[vertices] = new_node;
              }
          }
    };

      if (newly_refined_elements && !_mesh.is_replicated())
        {
          bool have_edge_queries = !edge_queries.empty();
          _mesh.comm().max(have_edge_queries);
          refinement_datum * refinement_data_ex = nullptr;
          if (have_edge_queries)
            Parallel::pull_parallel_vector_data
              (_mesh.comm(), edge_queries, edge_gather_functor,
               edge_action_functor, refinement_data_ex);
        }
    }

      // and elements that still have temporary id and unique_id
      // values.  Give them permanent ones.
      std::unordered_map<processor_id_type, std::vector<dof_id_type>> elems_to_query;
      for (const auto & [coarse_id, added_elem_map] : added_elements)
        {
          if (added_elem_map.empty())
            continue;

          const processor_id_type pid =
            added_elem_map.begin()->second->processor_id();
          if (pid == _mesh.processor_id())
            continue;

          elems_to_query[pid].push_back(coarse_id);
        }

      // Return fine element data based on vertex_average()

        elem_refinement_datum;

      auto added_gather_functor =
        [this]
        (processor_id_type,
         const std::vector<dof_id_type> & coarse_elems,
         std::vector<elem_refinement_datum> & coarse_refinements)
        {
          // Fill those requests
          const std::size_t query_size = coarse_elems.size();
          coarse_refinements.resize(query_size);
          for (auto i : make_range(query_size))
            {
              const dof_id_type coarse_id = coarse_elems[i];
              const auto & added =
                libmesh_map_find(added_elements, coarse_id);

              for (auto [vertex_avg, elem] : added)
                {
                  coarse_refinements[i].emplace_back
                    (vertex_avg,
                     elem->id()
#ifdef LIBMESH_ENABLE_UNIQUE_ID
                     , elem->unique_id()
#endif
                    );
                }
            }
        };

      auto added_action_functor =
        [this]
        (processor_id_type,
         const std::vector<dof_id_type> & coarse_elems,
         const std::vector<elem_refinement_datum> & coarse_refinements)
        {
          const std::size_t query_size = coarse_elems.size();
          for (auto i : make_range(query_size))
            {
              const dof_id_type coarse_id = coarse_elems[i];
              const auto & refinement_data = coarse_refinements[i];
              const auto & our_added =
                libmesh_map_find(added_elements, coarse_id);
              for (auto [vertex_avg, id
#ifdef LIBMESH_ENABLE_UNIQUE_ID
                         , unique_id
#endif
                   ] : refinement_data)
                {
                  Elem & our_elem = *libmesh_map_find(our_added,
                                                      vertex_avg);
                  libmesh_assert_equal_to(our_elem.vertex_average(),
                                          vertex_avg);

#ifdef LIBMESH_ENABLE_UNIQUE_ID
                  our_elem.set_unique_id(unique_id);
#endif
                  _mesh.renumber_elem(our_elem.id(), id);
                }
            }
        };

      elem_refinement_datum * refinement_data_ex = nullptr;
      Parallel::pull_parallel_vector_data
        (_mesh.comm(), elems_to_query, added_gather_functor,
         added_action_functor, refinement_data_ex);

      // That took care of our element ids; now get node ids.
      MeshCommunication mc;
      mc.make_node_proc_ids_parallel_consistent (_mesh);
      mc.make_node_ids_parallel_consistent (_mesh);
      mc.make_node_unique_ids_parallel_consistent (_mesh);
    }

  // In theory the operations on ghost elements are, after remote edge

}


void SimplexRefiner::fill_refinement_datum(std::pair<Node *, Node *> vertices,
                                           refinement_datum & vec)
{
  auto it = new_nodes.find(vertices);
  if (it == new_nodes.end())
    return;

  vec.emplace_back(vertices.first->id(),
                   vertices.second->id(),
                   it->second->id(),
                   it->second->processor_id());

  std::pair<Node *, Node *> subedge {vertices.first, it->second};
  if ((Point&)(*subedge.first) >
      (Point&)(*subedge.second))
    std::swap(subedge.first, subedge.second);
  fill_refinement_datum(subedge, vec);

  subedge = std::make_pair(it->second, vertices.second);
  if ((Point&)(*subedge.first) >
      (Point&)(*subedge.second))
    std::swap(subedge.first, subedge.second);
  fill_refinement_datum(subedge, vec);
}

          if (!mesh.is_replicated() &&
              hi_node->processor_id() != my_pid &&
              chosen_pid == my_pid)
            mesh.own_node(*hi_node);

          hi_node->processor_id() = chosen_pid;
        }

        // Hold off on trying to set the interior parent because we may actually
        // add lower dimensional elements before their interior parents
        if (old->interior_parent())
          ip_map[old] = el.get();

#ifdef LIBMESH_ENABLE_AMR
        if (old->has_children())

    // that aren't to that same third mesh.
    if (!ip_map.empty())
      {
        std::atomic<bool> existing_interior_parents{false};

        Threads::parallel_for
          (this->element_stored_range(),
           [&existing_interior_parents](const ElemRange & range)
           {
             for (Elem * elem : range)
               if (elem->interior_parent())
                 {
                   existing_interior_parents = true;
                   break;
                 }
           });

        MeshBase * other_interior_mesh =
          const_cast<MeshBase *>(&other_mesh.interior_mesh());

        // If we don't already have interior parents, then we can just
        // use whatever interior_mesh we need for the incoming
        // elements.
        if (!existing_interior_parents)
          {
            if (other_interior_mesh == &other_mesh)
              this->set_interior_mesh(*this);
            else
              this->set_interior_mesh(*other_interior_mesh);
          }

        if (other_interior_mesh == &other_mesh &&
            _interior_mesh == this)
          for (auto & elem_pair : ip_map)
            elem_pair.second->set_interior_parent(
              this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset));
        else if (other_interior_mesh == _interior_mesh)
          for (auto & elem_pair : ip_map)
            {

                       // did not, then those likewise must be remote
                       // children.
                       (current_elem->subactive() && neigh->has_children() &&
                        (neigh->level()+1) == current_elem->level())))
                    {
#ifdef DEBUG
                      // Let's make sure that "had children made remote"

                        libmesh_assert_not_equal_to (current_elem->processor_id(),
                                                     this->processor_id());
#endif // DEBUG
                      neigh = const_cast<RemoteElem *>(remote_elem);
                    }
                  // If neigh and current_elem are more than one level
                  // apart, figuring out whether we have a remote

          // children to search.
          if (pip == remote_elem || pip->active())
            {
              current_elem->set_interior_parent(pip);
              continue;
            }

          // For node comparisons we'll need a sensible tolerance

              if (&child == remote_elem)
                {
                  current_elem->set_interior_parent
                    (const_cast<RemoteElem *>(remote_elem));
                  continue;
                }

              bool child_contains_our_nodes = true;

      libmesh_fallthrough();
    case RIGHT_MULTIPLY_TRANSPOSE:
      {
        result_size = other.m() * this->m();
        if (other.n() == this->n())
          break;
      }
      libmesh_fallthrough();

  else if (flag == RIGHT_MULTIPLY_TRANSPOSE)
    {
      transa[0] = 't';
      std::swap(M,K);
    }

    case LEFT_MULTIPLY:            { this->_m = other.m(); break; }
    case RIGHT_MULTIPLY:           { this->_n = other.n(); break; }
    case LEFT_MULTIPLY_TRANSPOSE:  { this->_m = other.n(); break; }
    case RIGHT_MULTIPLY_TRANSPOSE: { this->_n = other.m(); break; }
    default:
      libmesh_error_msg("Unknown flag selected.");
    }

  // If we're not reusing submatrix and submatrix is already initialized
  // then we need to clear it, otherwise we get a memory leak.
  if (!reuse_submatrix && submatrix.initialized())
    submatrix.clear();

  // Construct row and column index sets.
  WrappedPetsc<IS> isrow;

T PetscMatrix<T>::operator () (const numeric_index_type i_in,
                               const numeric_index_type j_in) const
{
  libmesh_assert (this->initialized());

  // If the entry is not in the sparse matrix, it is 0.
  T value=0.;

  // this->close();
  // but closing it is a semiparallel operation; we want operator()
  // to run on one processor.
  libmesh_assert(this->closed());

  LibmeshPetscCall(MatGetValue(this->_mat, i_val, j_val, &value));

template <typename T>
void
PetscMatrixShellMatrix<T>::init(ParallelType libmesh_dbg_var(type))
{
#ifndef NDEBUG
  libmesh_assert(this->_dof_map);
  const auto m = this->_dof_map->n_dofs();
  const auto m_l = this->_dof_map->n_local_dofs();
  if (m != m_l)
    libmesh_assert(type != SERIAL);
#endif
  init_shell_mat(*this);
  this->set_context();
}

template <typename T>
void

    CasePCSetType(IDENTITY_PRECOND,     PCNONE)
    CasePCSetType(CHOLESKY_PRECOND,     PCCHOLESKY)
    CasePCSetType(ICC_PRECOND,          PCICC)
    CasePCSetType(ILU_PRECOND,          PCILU)
    CasePCSetType(LU_PRECOND,           PCLU)
    CasePCSetType(ASM_PRECOND,          PCASM)
    CasePCSetType(JACOBI_PRECOND,       PCJACOBI)

template <typename Obj>
void init_shell_mat(Obj & obj)
{
  libmesh_assert(obj._dof_map);

  numeric_index_type my_m = obj._dof_map->n_dofs();
  numeric_index_type m_l = obj._dof_map->n_local_dofs();
  if (obj._omit_constrained_dofs)
    {
      my_m -= obj._dof_map->n_constrained_dofs();
      m_l -= obj._dof_map->n_local_constrained_dofs();
    }

  const numeric_index_type my_n = my_m;
  const numeric_index_type n_l  = m_l;

  init_shell_mat(obj, my_m, my_n, m_l, n_l, obj._dof_map->block_size());
}

  //      (my_local_size == my_size))
  // But we do need to stay in sync for degenerate cases
  if (this->n_processors() == 1)
    return;


  // Build a parallel vector, initialize it with the local

  // only one processor
  if (n_processors() == 1)
    {
      v_local.resize(n);

      LibmeshPetscCall(VecGetArray (_vec, &values));

      for (PetscInt i=0; i<n; i++)
        v_local[i] = static_cast<Real>(values[i]);

      LibmeshPetscCall(VecRestoreArray (_vec, &values));
    }

  // otherwise multiple processors

{

template <typename T>
ReducedSparseMatrix<T>::ReducedSparseMatrix(const Parallel::Communicator & comm_in,
                                            const ReducedDofMap & reduced_dof_map)
  : PetscMatrixBase<T>(comm_in), _reduced_dof_map(reduced_dof_map)
{
  // The wrapped matrix owns the Mat
  this->_destroy_mat_on_exit = false;
}

template <typename T>
ReducedSparseMatrix<T>::~ReducedSparseMatrix() = default;

template <typename T>
void ReducedSparseMatrix<T>::adopt_reduced_matrix()
{
  this->_mat = _reduced->mat();
  this->_is_initialized = true;
  // Solvers recover this object from the Mat, so it, not the wrapped matrix, is the context
  this->set_context();
}

template <typename T>
void ReducedSparseMatrix<T>::init(ParallelType)
{
  libmesh_assert(_reduced_dof_map.initialized());

  const numeric_index_type n = _reduced_dof_map.n_dofs();
  const numeric_index_type n_local = _reduced_dof_map.n_local_dofs();

  _reduced = std::make_unique<PetscMatrix<T>>(this->comm());
  _reduced->init(
      n, n, n_local, n_local, _reduced_dof_map.reduced_n_nz(), _reduced_dof_map.reduced_n_oz());
  this->adopt_reduced_matrix();
}

template <typename T>
void ReducedSparseMatrix<T>::init(const numeric_index_type m,
                                  const numeric_index_type n,
                                  const numeric_index_type m_l,
                                  const numeric_index_type n_l,

                                  const numeric_index_type noz,
                                  const numeric_index_type blocksize)
{
  _reduced = std::make_unique<PetscMatrix<T>>(this->comm());
  _reduced->init(m, n, m_l, n_l, nnz, noz, blocksize);
  this->adopt_reduced_matrix();
}

template <typename T>
void ReducedSparseMatrix<T>::clear() noexcept
{
  _reduced.reset();
  this->_mat = nullptr;
  this->_is_initialized = false;
}

template <typename T>
PetscMatrix<T> & ReducedSparseMatrix<T>::reduced_matrix()
{
  libmesh_assert(_reduced);
  return *_reduced;
}

template <typename T>
SparseMatrix<T> & ReducedSparseMatrix<T>::operator=(const SparseMatrix<T> & v)
{
  const auto & other = cast_ref<const ReducedSparseMatrix<T> &>(v);
  libmesh_assert(&other._reduced_dof_map == &_reduced_dof_map);
  *_reduced = *other._reduced;
  return *this;
}

template <typename T>
void ReducedSparseMatrix<T>::to_reduced(const std::vector<numeric_index_type> & full,
                                        std::vector<PetscInt> & reduced) const
{
  reduced.resize(full.size());
  for (const auto i : index_range(full))
    {
      const dof_id_type r = _reduced_dof_map.find_reduced_dof(full[i]);
      reduced[i] = (r == DofObject::invalid_id) ? PetscInt(-1) : static_cast<PetscInt>(r);
    }
}

template <typename T>
void ReducedSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
                                        const std::vector<numeric_index_type> & rows,
                                        const std::vector<numeric_index_type> & cols)
{
  libmesh_assert(this->initialized());
  libmesh_assert_equal_to(rows.size(), dm.m());
  libmesh_assert_equal_to(cols.size(), dm.n());

  if (rows.empty() || cols.empty())
    return;

  this->to_reduced(rows, _rows);
  this->to_reduced(cols, _cols);
  LibmeshPetscCall(MatSetValues(this->_mat,
                                cast_int<PetscInt>(_rows.size()),
                                _rows.data(),
                                cast_int<PetscInt>(_cols.size()),

}

template <typename T>
void ReducedSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
                                        const std::vector<numeric_index_type> & dof_indices)
{
  this->add_matrix(dm, dof_indices, dof_indices);
}

template <typename T>
void ReducedSparseMatrix<T>::set(const numeric_index_type i,
                                 const numeric_index_type j,
                                 const T value)
{
  const dof_id_type ri = _reduced_dof_map.find_reduced_dof(i);
  const dof_id_type rj = _reduced_dof_map.find_reduced_dof(j);
  if (ri != DofObject::invalid_id && rj != DofObject::invalid_id)
    _reduced->set(ri, rj, value);
}

template <typename T>
void ReducedSparseMatrix<T>::add(const numeric_index_type i,
                                 const numeric_index_type j,
                                 const T value)
{
  const dof_id_type ri = _reduced_dof_map.find_reduced_dof(i);
  const dof_id_type rj = _reduced_dof_map.find_reduced_dof(j);
  if (ri != DofObject::invalid_id && rj != DofObject::invalid_id)
    _reduced->add(ri, rj, value);
}

template <typename T>
void ReducedSparseMatrix<T>::zero()
{
  _reduced->zero();
}

template <typename T>
void ReducedSparseMatrix<T>::close()
{
  _reduced->close();
}

template <typename T>
void ReducedSparseMatrix<T>::flush()
{
  _reduced->flush();
}

template <typename T>
std::unique_ptr<SparseMatrix<T>> ReducedSparseMatrix<T>::zero_clone() const
{
  auto ret = std::make_unique<ReducedSparseMatrix<T>>(this->comm(), _reduced_dof_map);
  ret->_reduced.reset(cast_ptr<PetscMatrix<T> *>(_reduced->zero_clone().release()));
  ret->adopt_reduced_matrix();
  return ret;
}

template <typename T>
std::unique_ptr<SparseMatrix<T>> ReducedSparseMatrix<T>::clone() const
{
  auto ret = std::make_unique<ReducedSparseMatrix<T>>(this->comm(), _reduced_dof_map);
  ret->_reduced.reset(cast_ptr<PetscMatrix<T> *>(_reduced->clone().release()));
  ret->adopt_reduced_matrix();
  return ret;
}

template <typename T>
void ReducedSparseMatrix<T>::add(const T a, const SparseMatrix<T> & X)
{
  const auto & other = cast_ref<const ReducedSparseMatrix<T> &>(X);
  _reduced->add(a, *other._reduced);
}

template <typename T>
T ReducedSparseMatrix<T>::operator()(const numeric_index_type i, const numeric_index_type j) const
{
  return (*_reduced)(i, j);
}

template <typename T>
Real ReducedSparseMatrix<T>::l1_norm() const
{
  return _reduced->l1_norm();
}

template <typename T>
Real ReducedSparseMatrix<T>::linfty_norm() const
{
  return _reduced->linfty_norm();
}

template <typename T>
void ReducedSparseMatrix<T>::print_personal(std::ostream & os) const
{
  _reduced->print_personal(os);
}

template <typename T>
void ReducedSparseMatrix<T>::get_diagonal(NumericVector<T> & dest) const
{
  _reduced->get_diagonal(dest);
}

template <typename T>
void ReducedSparseMatrix<T>::get_transpose(SparseMatrix<T> & dest) const
{
  auto & other = cast_ref<ReducedSparseMatrix<T> &>(dest);
  libmesh_assert(&other._reduced_dof_map == &_reduced_dof_map);
  if (!other._reduced)
    other._reduced = std::make_unique<PetscMatrix<T>>(this->comm());
  // The transpose replaces the destination's Mat, which the destination must then adopt
  _reduced->get_transpose(*other._reduced);
  other.adopt_reduced_matrix();
}

template <typename T>
void ReducedSparseMatrix<T>::get_row(numeric_index_type i,
                                     std::vector<numeric_index_type> & indices,
                                     std::vector<T> & values) const
{
  _reduced->get_row(i, indices, values);
}

template <typename T>
void ReducedSparseMatrix<T>::scale(const T scale)
{
  _reduced->scale(scale);
}

template class LIBMESH_EXPORT ReducedSparseMatrix<Number>;

  if (!this->initialized())
    {
      // The operator acts on the parent space of the static condensation dof map
      const auto n = _reduced_dof_map.parent_n_dofs();
      const auto n_local = _reduced_dof_map.parent_n_local_dofs();
      libmesh_assert(n == n_local || type != SERIAL);
      PetscMatrixShellMatrix<Number>::init(n, n, n_local, n_local);
      _parallel_type = type;
      this->init();
    }

  _reduced_rhs->init(n, n_local);

  // Initialize the reduced system matrix
  _sp = &_reduced_dof_map.reduced_sparsity();
  _reduced_sys_mat = SparseMatrix<Number>::build(this->comm());
  if (auto * const petsc_mat = dynamic_cast<PetscMatrix<Number> *>(_reduced_sys_mat.get()))
    {
      // Optimization for PETSc. This is critical for problems in which there are SCALAR dofs that
      // introduce dense rows to avoid allocating a dense matrix
      petsc_mat->init(
          n, n, n_local, n_local, _reduced_dof_map.reduced_n_nz(), _reduced_dof_map.reduced_n_oz());
    }
  else
    {

  // Build the ghosted parent-space solution vector. Note that this is, in general, *not equal* to
  // the system solution, e.g. this may correspond to the solution for the Newton *update*
  _ghosted_parent_sol = NumericVector<Number>::build(this->comm());
  _ghosted_parent_sol->init(_reduced_dof_map.parent_n_dofs(),
                            _reduced_dof_map.parent_n_local_dofs(),
                            _reduced_dof_map.parent_send_list(),
                            /*fast=*/false,
                            GHOSTED);

  _sc_is_initialized = true;
}

      for (const auto i : make_range(result.rows()))
        for (const auto j : make_range(result.cols()))
          shim(i, j) = result(i, j);
      for (const auto & var_reduced_space_indices : _reduced_dof_map.elem_dof_indices(elem_id))
        reduced_space_indices.insert(reduced_space_indices.end(),
                                     var_reduced_space_indices.begin(),
                                     var_reduced_space_indices.end());

void StaticCondensation::setup() { libmesh_assert(this->closed()); }

numeric_index_type StaticCondensation::m() const { return _reduced_dof_map.parent_n_dofs(); }

numeric_index_type StaticCondensation::row_start() const
{
  return _reduced_dof_map.parent_first_dof();
}

numeric_index_type StaticCondensation::row_stop() const
{
  return _reduced_dof_map.parent_end_dof();
}

void StaticCondensation::set(const numeric_index_type full_i,
                             const numeric_index_type full_j,
                             const Number val)
{
  const auto reduced_i = _reduced_dof_map.find_reduced_dof(full_i);
  const auto reduced_j = _reduced_dof_map.find_reduced_dof(full_j);
  libmesh_assert_msg((reduced_i != DofObject::invalid_id && reduced_j != DofObject::invalid_id) ||
                         _reduced_dof_map.parent(),
                     "Only uncondensed entries may be set in a statically condensed matrix");
  if (reduced_i != DofObject::invalid_id && reduced_j != DofObject::invalid_id)
    _reduced_sys_mat->set(reduced_i, reduced_j, val);
}

void StaticCondensation::set_current_elem(const Elem & elem)

  // Entries in rows or columns the parent space eliminates are dropped; the element matrix already
  // carries the parent reduction, e.g. C^T K C for eliminated constraints
  const auto * const parent = _reduced_dof_map.parent();
  auto in_parent = [parent](const auto global_index) {
    return !parent || parent->find_reduced_dof(global_index) != DofObject::invalid_id;
  };

  auto info_from_index = [&dof_data](const auto global_index) {
    auto index_it = dof_data.condensed_global_to_local_map.find(global_index);

      {
        const auto global_i = rows[i];
        const auto global_j = cols[j];
        if (!in_parent(global_i) || !in_parent(global_j))
          continue;
        const auto [i_is_condensed, local_i] = info_from_index(global_i);
        const auto [j_is_condensed, local_j] = info_from_index(global_j);
        if (i_is_condensed)

    elem_dof_values(i) = elem_dof_values_vec[i];
}

void StaticCondensation::forward_elimination(const NumericVector<Number> & parent_rhs)
{
  std::vector<Number> elem_condensed_rhs_vec;
  EigenVector elem_condensed_rhs, elem_uncondensed_rhs;

  parent_rhs.create_subvector(
      *_reduced_rhs, _reduced_dof_map.local_retained_parent_dofs(), /*all_global_entries=*/false);

  std::vector<dof_id_type> reduced_space_indices;
  for (auto elem : _mesh.active_local_element_ptr_range())

      auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
      reduced_space_indices.clear();
      const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());
      for (const auto & var_reduced_space_indices : _reduced_dof_map.elem_dof_indices(elem->id()))
        reduced_space_indices.insert(reduced_space_indices.end(),
                                     var_reduced_space_indices.begin(),
                                     var_reduced_space_indices.end());

      set_local_vectors(
          parent_rhs, dof_data.condensed_parent_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
      elem_uncondensed_rhs = -matrix_data.Auc * matrix_data.AccFactor.solve(elem_condensed_rhs);

      libmesh_assert(cast_int<std::size_t>(elem_uncondensed_rhs.size()) ==

  _reduced_rhs->close();
}

void StaticCondensation::backwards_substitution(const NumericVector<Number> & parent_rhs,
                                                NumericVector<Number> & parent_sol)
{
  std::vector<Number> elem_condensed_rhs_vec, elem_uncondensed_sol_vec;

      auto & matrix_data = libmesh_map_find(_elem_to_matrix_data, elem->id());
      const auto & dof_data = libmesh_map_find(_reduced_dof_map._elem_to_dof_data, elem->id());

      set_local_vectors(
          parent_rhs, dof_data.condensed_parent_dofs, elem_condensed_rhs_vec, elem_condensed_rhs);
      set_local_vectors(*_ghosted_parent_sol,
                        dof_data.uncondensed_parent_dofs,
                        elem_uncondensed_sol_vec,
                        elem_uncondensed_sol);

      elem_condensed_sol =
          matrix_data.AccFactor.solve(elem_condensed_rhs - matrix_data.Acu * elem_uncondensed_sol);
      parent_sol.insert(elem_condensed_sol.data(), dof_data.condensed_parent_dofs);
    }

  parent_sol.close();
}

void StaticCondensation::apply(const NumericVector<Number> & parent_rhs,
                               NumericVector<Number> & parent_sol)
{
  forward_elimination(parent_rhs);
  // Apparently PETSc will send us the yvec without zeroing it ahead of time. This can be a poor
  // initial guess for the Krylov solve as well as lead to bewildered users who expect their initial
  // residual norm to equal the norm of the RHS
  parent_sol.zero();
  _reduced_sol = parent_sol.get_subvector(_reduced_dof_map.local_retained_parent_dofs());
  _reduced_solver->solve(*_reduced_sys_mat, *_reduced_sol, *_reduced_rhs, 1e-5, 300);
  // Must restore to the parent-space solution because during backwards substitution we will need
  // to be able to read ghosted dofs and we don't support ghosting of subvectors
  parent_sol.restore_subvector(std::move(_reduced_sol),
                               _reduced_dof_map.local_retained_parent_dofs());
  *_ghosted_parent_sol = parent_sol;
  backwards_substitution(parent_rhs, parent_sol);
}

SolverPackage StaticCondensation::solver_package() { return libMesh::default_solver_package(); }

StaticCondensationDofMap::StaticCondensationDofMap(MeshBase & mesh,
                                                   System & system,
                                                   const DofMap & dof_map)
  : ReducedDofMap(mesh, system, dof_map, "reduced")
{
  libmesh_experimental();
}

StaticCondensationDofMap::~StaticCondensationDofMap() = default;

bool StaticCondensationDofMap::retain_field_dof(const Elem & elem,
                                                const unsigned int node_num,
                                                const unsigned int var_num,
                                                dof_id_type /*field_dof*/) const
{
  if (_uncondensed_vars.count(var_num))
    {
      libmesh_assert_msg(
          node_num == invalid_uint,
          "Users should not be providing continuous FEM variables to the uncondensed vars API");
      return true;
    }

  return node_num != invalid_uint && !elem.is_internal(node_num);
}

void StaticCondensationDofMap::retained_dof_added(const Elem & elem, const dof_id_type full_dof)
{
  auto & dof_data = _elem_to_dof_data[elem.id()];
  auto & map = dof_data.uncondensed_global_to_local_map;
  const auto local_number = cast_int<dof_id_type>(map.size());
  map.emplace(full_dof, local_number);
  dof_data.uncondensed_parent_dofs.push_back(this->parent_dof(full_dof));
}

void StaticCondensationDofMap::unretained_dof_found(const Elem & elem, const dof_id_type full_dof)
{
  auto & dof_data = _elem_to_dof_data[elem.id()];
  auto & map = dof_data.condensed_global_to_local_map;
  const auto local_number = cast_int<dof_id_type>(map.size());
  map.emplace(full_dof, local_number);
  dof_data.condensed_parent_dofs.push_back(this->parent_dof(full_dof));
}

void StaticCondensationDofMap::reinit()
{
  ReducedDofMap::reinit();

  // StaticCondensation expects an entry for every active local element, including those without
  // any degrees of freedom
  for (const auto * const elem : _mesh.active_local_element_ptr_range())
    _elem_to_dof_data[elem->id()];
}

void StaticCondensationDofMap::clear()
{
  ReducedDofMap::clear();
  _elem_to_dof_data.clear();
  _uncondensed_vars.clear();
}

template <>
void
WrappedPetsc<MatNullSpace>::destroy()
{
  LibmeshPetscCallA(PETSC_COMM_SELF, MatNullSpaceDestroy(&obj));
}

template <>
void

              mesh->interior_mesh().query_elem_ptr(interior_parent_id);
            if (!ip )
              elem->set_interior_parent
                (const_cast<RemoteElem *>(remote_elem));
            else
              elem->set_interior_parent(ip);
          }

              if ((e/blksize) < n)
                elem->processor_id() = cast_int<processor_id_type>(e/blksize);
              else
                elem->processor_id() = 0;
            }

          e++;

      // Otherwise  use kway
      else
        Metis::METIS_PartGraphKway(&n,
                                   &ncon,
                                   csr_graph.offsets.data(),
                                   csr_graph.vals.data(),

         // nodes.
         if (elem->type() == NODEELEM &&
             elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
           _pmetis->vwgt[local_index] = 50;
         else
           _pmetis->vwgt[local_index] = elem->n_nodes();

            data[i] = node.processor_id();
          }
        else
          data[i] = DofObject::invalid_processor_id;
      }
  }

                bad_pids.erase(it);
              }
            else
              proc_id = node.choose_processor_id(proc_id, new_proc_id);

            if (proc_id == new_proc_id)
              data_changed = true;

          // Use a set to avoid duplicates.  Use a well-defined method
          // of ordering that set to make debugging easier.
          std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
          for (const Node & node : elem->node_ref_range())
            {
              if (const auto row_it = mesh_constrained_nodes.find(&node);
                  row_it != end_it)
                for (const auto & [pr, coef] : row_it->second)
                  {
                    libmesh_ignore(coef); // avoid gcc 7 warning
                    constraining_elems.insert(pr.first);
                  }
            }
          for (const Elem * constraining_elem : constraining_elems)
            elems_constrained_by.emplace(constraining_elem, elem);
        }

        {
          // get all relevant interior elements
          std::set<const Elem *> neighbor_set;
          elem->find_interior_neighbors(neighbor_set);

          for (const auto & neighbor : neighbor_set)
            interior_to_boundary_map.emplace(neighbor, elem);
        }
    }

        {
          // get all relevant interior elements
          std::set<const Elem *> neighbor_set;
          elem->find_interior_neighbors(neighbor_set);

          for (const auto & neighbor : neighbor_set)
            {
              const dof_id_type neighbor_global_index_by_pid =
                _global_index_by_pid_map[neighbor->id()];

              graph_row.push_back(neighbor_global_index_by_pid);
            }
        }

      // Check for any boundary neighbors
      for (const auto & pr : as_range(interior_to_boundary_map.equal_range(elem)))
        {
          const Elem * neighbor = pr.second;

          const dof_id_type neighbor_global_index_by_pid =
            _global_index_by_pid_map[neighbor->id()];

          graph_row.push_back(neighbor_global_index_by_pid);
        }

      // Check for any constraining elements

          // Use a set to avoid duplicates.  Use a well-defined method
          // of ordering that set to make debugging easier.
          std::set<const Elem *, CompareElemIdsByLevel> constraining_elems;
          for (const Node & node : elem->node_ref_range())
            {
              if (const auto row_it = mesh_constrained_nodes.find(&node);
                  row_it != end_it)
                for (const auto & [pr, coef] : row_it->second)
                  {
                    libmesh_ignore(coef); // avoid gcc 7 warning
                    constraining_elems.insert(pr.first);
                  }
            }
          for (const Elem * constraining_elem : constraining_elems)
            {
              const dof_id_type constraining_global_index_by_pid =
                _global_index_by_pid_map[constraining_elem->id()];

              graph_row.push_back(constraining_global_index_by_pid);

              // We can't be sure if the constraining element's owner sees
              // the assembly element, so to get a symmetric connectivity
              // graph we'll need to tell them about us to be safe.
              if (constraining_elem->processor_id() != mesh.processor_id())
                connections_to_push[constraining_elem->processor_id()].emplace_back
                  (global_index_by_pid, constraining_global_index_by_pid);
            }
        }

      // Check for any constrained elements
      for (const auto & pr : as_range(elems_constrained_by.equal_range(elem)))
        {
          const Elem * constrained = pr.second;
          const dof_id_type constrained_global_index_by_pid =
            _global_index_by_pid_map[constrained->id()];

          graph_row.push_back(constrained_global_index_by_pid);

          // We can't be sure if the constrained element's owner sees
          // the assembly element, so to get a symmetric connectivity
          // graph we'll need to tell them about us to be safe.
          if (constrained->processor_id() != mesh.processor_id())
            connections_to_push[constrained->processor_id()].emplace_back
              (global_index_by_pid, constrained_global_index_by_pid);
        }
    }

  // processor doesn't see our assembly element.  Let's push those
  // entries, to ensure they're counted on both sides.
  auto symmetrize_entries =
    [this, first_local_elem]
    (processor_id_type /*src_pid*/,
     const std::vector<std::pair<dof_id_type, dof_id_type>> & incoming_entries)
    {
      for (auto [i, j] : incoming_entries)
        {
          libmesh_assert_greater_equal(j, first_local_elem);
          const std::size_t jl = j - first_local_elem;
          libmesh_assert_less(jl, _dual_graph.size());
          std::vector<dof_id_type> & graph_row = _dual_graph[jl];
          if (std::find(graph_row.begin(), graph_row.end(), i) == graph_row.end())
          {
//            std::cerr << "Pushing back (" << j << ", " << i << ") from " << src_pid << std::endl;
            graph_row.push_back(i);
          }
        }
    };

      // periodic constraints
      if (assemble_matrix && symmetrize)
        {
          DenseMatrix<Number> Ke_transpose;
          context.get_elem_jacobian().get_transpose(Ke_transpose);
          context.get_elem_jacobian() += Ke_transpose;
          context.get_elem_jacobian() *= 0.5;
        }

      // As discussed above, we can set apply_dof_constraints=false to
      // get A instead of C^T A C

                               apply_dof_constraints);
}

void RBConstruction::add_scaled_Aq(Number scalar,
                                   unsigned int q_a,
                                   SparseMatrix<Number> * input_matrix,
                                   bool symmetrize)
{
  LOG_SCOPE("add_scaled_Aq()", "RBConstruction");

  libmesh_error_msg_if(q_a >= get_rb_theta_expansion().get_n_A_terms(),
                       "Error: We must have q < Q_a in add_scaled_Aq.");

  if (!symmetrize)
    {
      input_matrix->add(scalar, *get_Aq(q_a));
      input_matrix->close();
    }
  else
    {
      add_scaled_matrix_and_vector(scalar,
                                   &rb_assembly_expansion->get_A_assembly(q_a),
                                   input_matrix,
                                   nullptr,
                                   symmetrize);
    }
}

void RBConstruction::assemble_misc_matrices()
{

  if (rel_greedy_error < this->rel_training_tolerance)
    {
      libMesh::out << "Relative error tolerance reached." << std::endl;
      return true;
    }

  RBEvaluation & rbe = get_rb_evaluation();

  unsigned int remainder = n_global_samples%communicator.size();
  if (communicator.rank() < remainder)
    {
      n_local_samples = (quotient + 1);
      first_local_index = communicator.rank()*(quotient+1);
    }
  else
    {

}

template <class Base>
void RBConstructionBase<Base>::set_params_from_training_set_and_broadcast(unsigned int global_index)
{
  libmesh_error_msg_if(!_training_parameters_initialized,
                       "Error: training parameters must first be initialized.");

  processor_id_type root_id = 0;
  if ((this->get_first_local_training_index() <= global_index) &&
      (global_index < this->get_last_local_training_index()))
    {
      // Set parameters on only one processor
      set_params_from_training_set(global_index);

      // set root_id, only non-zero on one processor
      root_id = this->processor_id();
    }

  // broadcast
  this->comm().max(root_id);
  broadcast_parameters(root_id);
}

template <class Base>
void RBConstructionBase<Base>::initialize_training_parameters(const RBParameters & mu_min,

    }
  else
    {
      if (!serial_training_set)
        {
          // seed the random number generator with the provided value
          // and the processor ID so that the seed is different
          // on different processors
          std::srand( static_cast<unsigned>( training_parameters_random_seed*(1+communicator.rank()) ));
        }
      else
        {

namespace libMesh
{

RBSCMConstruction::RBSCMConstruction (EquationSystems & es,
                                      const std::string & name_in,
                                      const unsigned int number_in)
  : Parent(es, name_in, number_in),
    SCM_training_tolerance(0.5),
    RB_system_name(""),
    rb_scm_eval(nullptr)
{

  // set assemble_before_solve flag to false
  // so that we control matrix assembly.
  assemble_before_solve = false;

  // We symmetrize all operators hence use symmetric solvers
  set_eigenproblem_type(GHEP);

}

RBSCMConstruction::~RBSCMConstruction () = default;

void RBSCMConstruction::clear()
{
  Parent::clear();
}

void RBSCMConstruction::set_rb_scm_evaluation(RBSCMEvaluation & rb_scm_eval_in)
{
  rb_scm_eval = &rb_scm_eval_in;
}

RBSCMEvaluation & RBSCMConstruction::get_rb_scm_evaluation()
{
  libmesh_error_msg_if(!rb_scm_eval, "Error: RBSCMEvaluation object hasn't been initialized yet");

  return *rb_scm_eval;
}

RBThetaExpansion & RBSCMConstruction::get_rb_theta_expansion()
{
  return get_rb_scm_evaluation().get_rb_theta_expansion();
}

void RBSCMConstruction::process_parameters_file(const std::string & parameters_filename)
{
  // First read in data from parameters_filename
  GetPot infile(parameters_filename);
  const unsigned int n_training_samples = infile("n_training_samples",1);
  const bool deterministic_training     = infile("deterministic_training",false);

  // Read in training_parameters_random_seed value.  This is used to
  // seed the RNG when picking the training parameters.  By default the
  // value is -1, which means use std::time to seed the RNG.
  unsigned int training_parameters_random_seed_in = static_cast<unsigned int>(-1);
  training_parameters_random_seed_in = infile("training_parameters_random_seed",
                                              training_parameters_random_seed_in);
  set_training_random_seed(static_cast<int>(training_parameters_random_seed_in));

  // SCM Greedy termination tolerance
  const Real SCM_training_tolerance_in = infile("SCM_training_tolerance", SCM_training_tolerance);
  set_SCM_training_tolerance(SCM_training_tolerance_in);

  // Initialize the parameter ranges and the parameters themselves
  unsigned int n_continuous_parameters = infile.vector_variable_size("parameter_names");
  RBParameters mu_min_in;
  RBParameters mu_max_in;
  for (unsigned int i=0; i<n_continuous_parameters; i++)
    {
      // Read in the parameter names
      std::string param_name = infile("parameter_names", "NONE", i);

      {
        Real min_val = infile(param_name, 0., 0);
        mu_min_in.set_value(param_name, min_val);
      }

      {
        Real max_val = infile(param_name, 0., 1);
        mu_max_in.set_value(param_name, max_val);
      }
    }

  std::map<std::string, std::vector<Real>> discrete_parameter_values_in;

  unsigned int n_discrete_parameters = infile.vector_variable_size("discrete_parameter_names");
  for (unsigned int i=0; i<n_discrete_parameters; i++)
    {
      std::string param_name = infile("discrete_parameter_names", "NONE", i);

      discrete_parameter_values_in[param_name] = vals_for_param;
    }

  initialize_parameters(mu_min_in, mu_max_in, discrete_parameter_values_in);

  std::map<std::string,bool> log_scaling;
  const RBParameters & mu = get_parameters();
  unsigned int i=0;
  for (const auto & pr : mu)
    {
      const std::string & param_name = pr.first;
      log_scaling[param_name] = static_cast<bool>(infile("log_scaling", 0, i++));
    }

  initialize_training_parameters(this->get_parameters_min(),
                                 this->get_parameters_max(),
                                 n_training_samples,
                                 log_scaling,
                                 deterministic_training);   // use deterministic parameters
}

void RBSCMConstruction::print_info()
{
  // Print out info that describes the current setup
  libMesh::out << std::endl << "RBSCMConstruction parameters:" << std::endl;
  libMesh::out << "system name: " << this->name() << std::endl;
  libMesh::out << "SCM Greedy tolerance: " << get_SCM_training_tolerance() << std::endl;
  if (rb_scm_eval)
    {
      libMesh::out << "A_q operators attached: " << get_rb_theta_expansion().get_n_A_terms() << std::endl;
    }
  else
    {
      libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
    }
  libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
  for (const auto & pr : get_parameters())
    {
      const std::string & param_name = pr.first;
      libMesh::out <<   "Parameter " << param_name
                   << ": Min = " << get_parameter_min(param_name)
                   << ", Max = " << get_parameter_max(param_name) << std::endl;
    }
  print_discrete_parameter_values();
  libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
  libMesh::out << std::endl;
}

void RBSCMConstruction::resize_SCM_vectors()
{

  rb_scm_eval->B_max.resize(get_rb_theta_expansion().get_n_A_terms());
}

void RBSCMConstruction::add_scaled_symm_Aq(unsigned int q_a, Number scalar)
{
  LOG_SCOPE("add_scaled_symm_Aq()", "RBSCMConstruction");
  // Load the operators from the RBConstruction
  EquationSystems & es = this->get_equation_systems();
  RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);
  rb_system.add_scaled_Aq(scalar, q_a, &get_matrix_A(), true);
}

void RBSCMConstruction::load_matrix_B()
{
  // Load the operators from the RBConstruction
  EquationSystems & es = this->get_equation_systems();
  RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);

  matrix_B->zero();
  matrix_B->close();
  matrix_B->add(1.,*rb_system.get_inner_product_matrix());
}

void RBSCMConstruction::perform_SCM_greedy()
{
  LOG_SCOPE("perform_SCM_greedy()", "RBSCMConstruction");

  // initialize rb_scm_eval's parameters
  rb_scm_eval->initialize_parameters(*this);

#ifdef LIBMESH_ENABLE_CONSTRAINTS
  // Get a list of constrained dofs from rb_system
  std::set<dof_id_type> constrained_dofs_set;
  EquationSystems & es = this->get_equation_systems();
  RBConstruction & rb_system = es.get_system<RBConstruction>(RB_system_name);

  for (dof_id_type i=0; i<rb_system.n_dofs(); i++)
    {
      if (rb_system.get_dof_map().is_constrained_dof(i))
        {
          constrained_dofs_set.insert(i);
        }
    }

  // Use these constrained dofs to identify which dofs we want to "get rid of"
  // (i.e. condense) in our eigenproblems.
  this->initialize_condensed_dofs(constrained_dofs_set);
#endif // LIBMESH_ENABLE_CONSTRAINTS

  // Copy the inner product matrix over from rb_system to be used as matrix_B
  load_matrix_B();

  attach_deflation_space();

  compute_SCM_bounding_box();
  // This loads the new parameter into current_parameters
  enrich_C_J(0);

  unsigned int SCM_iter=0;
  while (true)
    {
      // matrix_A is reinitialized for the current parameters
      // on each call to evaluate_stability_constant
      evaluate_stability_constant();

      std::pair<unsigned int,Real> SCM_error_pair = compute_SCM_bounds_on_training_set();

      libMesh::out << "SCM iteration " << SCM_iter
                   << ", max_SCM_error = " << SCM_error_pair.second << std::endl;

      if (SCM_error_pair.second < SCM_training_tolerance)
        {
          libMesh::out << std::endl << "SCM tolerance of " << SCM_training_tolerance << " reached."
                       << std::endl << std::endl;
          break;
        }

      // If we need another SCM iteration, then enrich C_J
      enrich_C_J(SCM_error_pair.first);

      libMesh::out << std::endl << "-----------------------------------" << std::endl << std::endl;

      SCM_iter++;
    }
}

void RBSCMConstruction::compute_SCM_bounding_box()
{
  LOG_SCOPE("compute_SCM_bounding_box()", "RBSCMConstruction");

  // Resize the bounding box vectors
  rb_scm_eval->B_min.resize(get_rb_theta_expansion().get_n_A_terms());
  rb_scm_eval->B_max.resize(get_rb_theta_expansion().get_n_A_terms());

  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
    {
      matrix_A->zero();
      add_scaled_symm_Aq(q, 1.);

      // Compute B_min(q)
      eigen_solver->set_position_of_spectrum(SMALLEST_REAL);
      set_eigensolver_properties(q);

      solve();
      // TODO: Assert convergence for eigensolver

      unsigned int nconv = get_n_converged();
      if (nconv != 0)
        {
          std::pair<Real, Real> eval = get_eigenpair(0);

          // ensure that the eigenvalue is real
          libmesh_assert_less (eval.second, TOLERANCE);

          rb_scm_eval->set_B_min(q, eval.first);
          libMesh::out << std::endl << "B_min("<<q<<") = " << rb_scm_eval->get_B_min(q) << std::endl;
        }
      else
        libmesh_error_msg("Eigen solver for computing B_min did not converge");

      // Compute B_max(q)
      eigen_solver->set_position_of_spectrum(LARGEST_REAL);
      set_eigensolver_properties(q);

      solve();
      // TODO: Assert convergence for eigensolver

      nconv = get_n_converged();
      if (nconv != 0)
        {
          std::pair<Real, Real> eval = get_eigenpair(0);

          // ensure that the eigenvalue is real
          libmesh_assert_less (eval.second, TOLERANCE);

          rb_scm_eval->set_B_max(q,eval.first);
          libMesh::out << "B_max("<<q<<") = " << rb_scm_eval->get_B_max(q) << std::endl;
        }
      else
        libmesh_error_msg("Eigen solver for computing B_max did not converge");
    }
}

void RBSCMConstruction::evaluate_stability_constant()
{
  LOG_SCOPE("evaluate_stability_constant()", "RBSCMConstruction");

  // Get current index of C_J
  const unsigned int j = cast_int<unsigned int>(rb_scm_eval->C_J.size()-1);

  eigen_solver->set_position_of_spectrum(SMALLEST_REAL);

  // For non-coercive time-dependent problems, B is set to the mass matrix

  // Set matrix A corresponding to mu_star
  matrix_A->zero();
  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
    {
      add_scaled_symm_Aq(q, get_rb_theta_expansion().eval_A_theta(q,get_parameters()));
    }

  set_eigensolver_properties(-1);
  solve();
  // TODO: Assert convergence for eigensolver

  unsigned int nconv = get_n_converged();
  if (nconv != 0)
    {
      std::pair<Real, Real> eval = get_eigenpair(0);

      // ensure that the eigenvalue is real
      libmesh_assert_less (eval.second, TOLERANCE);

      // Store the coercivity constant corresponding to mu_star
      rb_scm_eval->set_C_J_stability_constraint(j,eval.first);
      libMesh::out << std::endl << "Stability constant for C_J("<<j<<") = "
                   << rb_scm_eval->get_C_J_stability_constraint(j) << std::endl << std::endl;

      // Compute and store the vector y = (y_1, \ldots, y_Q) for the
      // eigenvector currently stored in eigen_system.solution.
      // We use this later to compute the SCM upper bounds.
      Real norm_B2 = libmesh_real( B_inner_product(*solution, *solution) );

      for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
        {
          Real norm_Aq2 = libmesh_real( Aq_inner_product(q, *solution, *solution) );

          rb_scm_eval->set_SCM_UB_vector(j,q,norm_Aq2/norm_B2);
        }
    }
  else
    libmesh_error_msg("Error: Eigensolver did not converge in evaluate_stability_constant");
}

Number RBSCMConstruction::B_inner_product(const NumericVector<Number> & v,
                                          const NumericVector<Number> & w) const
{
  matrix_B->vector_mult(*inner_product_storage_vector, w);

  return v.dot(*inner_product_storage_vector);
}

Number RBSCMConstruction::Aq_inner_product(unsigned int q,
                                           const NumericVector<Number> & v,
                                           const NumericVector<Number> & w)
{
  libmesh_error_msg_if(q >= get_rb_theta_expansion().get_n_A_terms(),
                       "Error: We must have q < Q_a in Aq_inner_product.");

  matrix_A->zero();
  add_scaled_symm_Aq(q, 1.);
  matrix_A->vector_mult(*inner_product_storage_vector, w);

  return v.dot(*inner_product_storage_vector);
}

std::pair<unsigned int,Real> RBSCMConstruction::compute_SCM_bounds_on_training_set()
{
  LOG_SCOPE("compute_SCM_bounds_on_training_set()", "RBSCMConstruction");

  unsigned int new_C_J_index = 0;
  Real max_SCM_error = 0.;

  numeric_index_type first_index = get_first_local_training_index();
  for (unsigned int i=0; i<get_local_n_training_samples(); i++)
    {
      set_params_from_training_set(first_index+i);
      rb_scm_eval->set_parameters( get_parameters() );
      Real LB = rb_scm_eval->get_SCM_LB();
      Real UB = rb_scm_eval->get_SCM_UB();

      Real error_i = SCM_greedy_error_indicator(LB, UB);

      if (error_i > max_SCM_error)
        {
          max_SCM_error = error_i;
          new_C_J_index = i;
        }
    }

  numeric_index_type global_index = first_index + new_C_J_index;
  std::pair<numeric_index_type,Real> error_pair(global_index, max_SCM_error);
  get_global_max_error_pair(this->comm(),error_pair);

  return error_pair;
}

void RBSCMConstruction::enrich_C_J(unsigned int new_C_J_index)
{
  LOG_SCOPE("enrich_C_J()", "RBSCMConstruction");

  set_params_from_training_set_and_broadcast(new_C_J_index);

  rb_scm_eval->C_J.push_back(get_parameters());

  libMesh::out << std::endl << "SCM: Added mu = (";

  bool first = true;
  for (const auto & pr : get_parameters())
    {
      if (!first)
        libMesh::out << ",";
      const std::string & param_name = pr.first;
      RBParameters C_J_params = rb_scm_eval->C_J[rb_scm_eval->C_J.size()-1];
      libMesh::out << C_J_params.get_value(param_name);
      first = false;
    }
  libMesh::out << ")" << std::endl;

  // Finally, resize C_J_stability_vector and SCM_UB_vectors
  rb_scm_eval->C_J_stability_vector.push_back(0.);

  std::vector<Real> zero_vector(get_rb_theta_expansion().get_n_A_terms());
  rb_scm_eval->SCM_UB_vectors.push_back(zero_vector);
}


} // namespace libMesh

namespace libMesh
{

RBSCMEvaluation::RBSCMEvaluation (const Parallel::Communicator & comm_in) :
  ParallelObject(comm_in)
{
  // Clear SCM data vectors
  B_min.clear();

  C_J.clear();
  C_J_stability_vector.clear();
  SCM_UB_vectors.clear();
}

RBSCMEvaluation::~RBSCMEvaluation () = default;

void RBSCMEvaluation::set_rb_theta_expansion(RBThetaExpansion & rb_theta_expansion_in)
{
  rb_theta_expansion = &rb_theta_expansion_in;
}

RBThetaExpansion & RBSCMEvaluation::get_rb_theta_expansion()
{
  libmesh_error_msg_if(!rb_theta_expansion, "Error: rb_theta_expansion hasn't been initialized yet");

  return *rb_theta_expansion;
}

void RBSCMEvaluation::set_C_J_stability_constraint(unsigned int j, Real stability_const_in)
{
  libmesh_error_msg_if(j >= C_J_stability_vector.size(), "Error: Input parameter j is too large in set_C_J_stability_constraint.");

  // we assume that C_J_stability_vector is resized elsewhere
  // to be the same size as C_J.
  libmesh_assert_equal_to (C_J_stability_vector.size(), C_J.size());

  C_J_stability_vector[j] = stability_const_in;
}

Real RBSCMEvaluation::get_C_J_stability_constraint(unsigned int j) const
{
  libmesh_error_msg_if(j >= C_J_stability_vector.size(), "Error: Input parameter j is too large in get_C_J_stability_constraint.");

  return C_J_stability_vector[j];
}

void RBSCMEvaluation::set_SCM_UB_vector(unsigned int j, unsigned int q, Real y_q)
{
  // First make sure that j <= J
  libmesh_error_msg_if(j >= SCM_UB_vectors.size(), "Error: We must have j < J in set_SCM_UB_vector.");

  // Next make sure that q <= Q_a or Q_a_hat
  libmesh_error_msg_if(q >= SCM_UB_vectors[0].size(), "Error: q is too large in set_SCM_UB_vector.");

  SCM_UB_vectors[j][q] = y_q;
}

Real RBSCMEvaluation::get_SCM_UB_vector(unsigned int j, unsigned int q)
{
  // First make sure that j <= J
  libmesh_error_msg_if(j >= SCM_UB_vectors.size(), "Error: We must have j < J in get_SCM_UB_vector.");
  libmesh_error_msg_if(q >= SCM_UB_vectors[0].size(), "Error: q is too large in get_SCM_UB_vector.");

  return SCM_UB_vectors[j][q];
}

const RBParameters & RBSCMEvaluation::get_C_J_entry(unsigned int j)

  return C_J[j];
}

Real RBSCMEvaluation::get_B_min(unsigned int q) const
{
  libmesh_error_msg_if(q >= B_min.size(), "Error: q is too large in get_B_min.");

  return B_min[q];
}


Real RBSCMEvaluation::get_B_max(unsigned int q) const
{
  libmesh_error_msg_if(q >= B_max.size(), "Error: q is too large in get_B_max.");

  return B_max[q];
}

void RBSCMEvaluation::set_B_min(unsigned int q, Real B_min_val)
{
  libmesh_error_msg_if(q >= B_min.size(), "Error: q is too large in set_B_min.");

  B_min[q] = B_min_val;
}

void RBSCMEvaluation::set_B_max(unsigned int q, Real B_max_val)
{
  libmesh_error_msg_if(q >= B_max.size(), "Error: q is too large in set_B_max.");

  B_max[q] = B_max_val;
}

Real RBSCMEvaluation::get_SCM_LB()
{
  LOG_SCOPE("get_SCM_LB()", "RBSCMEvaluation");

  // Initialize the LP
  glp_prob * lp;
  lp = glp_create_prob();
  glp_set_obj_dir(lp,GLP_MIN);

  // Add columns to the LP: corresponds to
  // the variables y_1,...y_Q_a.
  // These are the same for each \mu in the SCM
  // training set, hence can do this up front.
  glp_add_cols(lp,rb_theta_expansion->get_n_A_terms());

  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
    {
      if (B_max[q] < B_min[q]) // Invalid bound, set as free variable
        {
          // GLPK indexing is not zero based!
          glp_set_col_bnds(lp, q+1, GLP_FR, 0., 0.);

      else
        {
          // GLPK indexing is not zero based!
          glp_set_col_bnds(lp, q+1, GLP_DB, double(B_min[q]), double(B_max[q]));
        }

      // If B_max is not defined, just set lower bounds...

  // variables that define the constraints at each
  // mu \in C_J_M
  unsigned int n_rows = cast_int<unsigned int>(C_J.size());
  glp_add_rows(lp, n_rows);

  // Now put current_parameters in saved_parameters
  save_current_parameters();

  unsigned int matrix_size = n_rows*rb_theta_expansion->get_n_A_terms();
  std::vector<int> ia(matrix_size+1);
  std::vector<int> ja(matrix_size+1);
  std::vector<double> ar(matrix_size+1);
  unsigned int count=0;
  for (unsigned int m=0; m<n_rows; m++)
    {
      set_current_parameters_from_C_J(m);

      // Set the lower bound on the auxiliary variable
      // due to the stability constant at mu_index
      glp_set_row_bnds(lp, m+1, GLP_LO, double(C_J_stability_vector[m]), 0.);

      // Now define the matrix that relates the y's
      // to the auxiliary variables at the current
      // value of mu.
      for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
        {
          count++;

          ia[count] = m+1;
          ja[count] = q+1;

          // This can only handle Reals right now
          ar[count] = double(libmesh_real( rb_theta_expansion->eval_A_theta(q,get_parameters()) ));
        }
    }

  // Now load the original parameters back into current_parameters
  // in order to set the coefficients of the objective function
  reload_current_parameters();

  glp_load_matrix(lp, matrix_size, ia.data(), ja.data(), ar.data());

  for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
    {
      glp_set_obj_coef(lp,q+1, double(libmesh_real( rb_theta_expansion->eval_A_theta(q,get_parameters()) )) );
    }

  // Use this command to initialize the basis for the LP

  //lpx_cpx_basis(lp); //glp_cpx_basis(lp);

  glp_smcp parm;
  glp_init_smcp(&parm);
  parm.msg_lev = GLP_MSG_ERR;
  parm.meth = GLP_DUAL;


  // use the simplex method and solve the LP
  glp_simplex(lp, &parm);

  Real min_J_obj = glp_get_obj_val(lp);

  //   int simplex_status =  glp_get_status(lp);
  //   if (simplex_status == GLP_UNBND)

  //   }

  // Destroy the LP
  glp_delete_prob(lp);

  return min_J_obj;
}

Real RBSCMEvaluation::get_SCM_UB()
{
  LOG_SCOPE("get_SCM_UB()", "RBSCMEvaluation");

  // to C_J_M (SCM_UB_vectors contains vectors for all of
  // C_J).
  Real min_J_obj = 0.;
  for (unsigned int m=0; m<n_rows; m++)
    {
      const std::vector<Real> UB_vector = SCM_UB_vectors[m];

      Real J_obj = 0.;
      for (unsigned int q=0; q<rb_theta_expansion->get_n_A_terms(); q++)
        {
          J_obj += libmesh_real( rb_theta_expansion->eval_A_theta(q,get_parameters()) )*UB_vector[q];
        }

      if ((m==0) || (J_obj < min_J_obj))
        {
          min_J_obj = J_obj;
        }
    }

  return min_J_obj;
}

void RBSCMEvaluation::set_current_parameters_from_C_J(unsigned int C_J_index)
{
  set_parameters(C_J[C_J_index]);
}

void RBSCMEvaluation::save_current_parameters()
{
  saved_parameters = get_parameters();
}

void RBSCMEvaluation::reload_current_parameters()
{
  set_parameters(saved_parameters);
}

void RBSCMEvaluation::legacy_write_offline_data_to_files(const std::string & directory_name,
                                                         const bool write_binary_data)
{
  LOG_SCOPE("legacy_write_offline_data_to_files()", "RBSCMEvaluation");

  if (this->processor_id() == 0)
    {
      // Make a directory to store all the data files
      if (mkdir(directory_name.c_str(), 0777) == -1)
        {
          libMesh::out << "In RBSCMEvaluation::write_offline_data_to_files, directory "
                       << directory_name << " already exists, overwriting contents." << std::endl;
        }

      // The writing mode: ENCODE for binary, WRITE for ASCII
      XdrMODE mode = write_binary_data ? ENCODE : WRITE;

      // The suffix to use for all the files that are written out
      const std::string suffix = write_binary_data ? ".xdr" : ".dat";

      // Stream for building the file names
      std::ostringstream file_name;

      // Write out the parameter ranges
      file_name.str("");
      file_name << directory_name << "/parameter_ranges" << suffix;
      std::string continuous_param_file_name = file_name.str();

      // Write out the discrete parameter values
      file_name.str("");
      file_name << directory_name << "/discrete_parameter_values" << suffix;
      std::string discrete_param_file_name = file_name.str();

      write_parameter_data_to_files(continuous_param_file_name,
                                    discrete_param_file_name,
                                    write_binary_data);

      // Write out the bounding box min values
      file_name.str("");
      file_name << directory_name << "/B_min" << suffix;
      Xdr B_min_out(file_name.str(), mode);

      for (auto i : make_range(B_min.size()))
        {
          Real B_min_i = get_B_min(i);
          B_min_out << B_min_i;
        }
      B_min_out.close();


      // Write out the bounding box max values
      file_name.str("");
      file_name << directory_name << "/B_max" << suffix;
      Xdr B_max_out(file_name.str(), mode);

      for (auto i : make_range(B_max.size()))
        {
          Real B_max_i = get_B_max(i);
          B_max_out << B_max_i;
        }
      B_max_out.close();

      // Write out the length of the C_J data
      file_name.str("");
      file_name << directory_name << "/C_J_length" << suffix;
      Xdr C_J_length_out(file_name.str(), mode);

      unsigned int C_J_length = cast_int<unsigned int>(C_J.size());
      C_J_length_out << C_J_length;
      C_J_length_out.close();

      // Write out C_J_stability_vector
      file_name.str("");
      file_name << directory_name << "/C_J_stability_vector" << suffix;
      Xdr C_J_stability_vector_out(file_name.str(), mode);

      for (auto i : make_range(C_J_stability_vector.size()))
        {
          Real C_J_stability_constraint_i = get_C_J_stability_constraint(i);
          C_J_stability_vector_out << C_J_stability_constraint_i;
        }
      C_J_stability_vector_out.close();

      // Write out C_J
      file_name.str("");
      file_name << directory_name << "/C_J" << suffix;
      Xdr C_J_out(file_name.str(), mode);

      for (const auto & param : C_J)
        for (const auto & pr : param)
          for (const auto & value_vector : pr.second)
            {
              // Need to make a copy of the value so that it's not const
              // Xdr is not templated on const's
              libmesh_error_msg_if(value_vector.size() != 1,
                                   "Error: multi-value RB parameters are not yet supported here.");
              Real param_value = value_vector[0];
              C_J_out << param_value;
            }
      C_J_out.close();

      // Write out SCM_UB_vectors get_SCM_UB_vector
      file_name.str("");
      file_name << directory_name << "/SCM_UB_vectors" << suffix;
      Xdr SCM_UB_vectors_out(file_name.str(), mode);

      for (auto i : make_range(SCM_UB_vectors.size()))
        for (auto j : make_range(rb_theta_expansion->get_n_A_terms()))
          {
            Real SCM_UB_vector_ij = get_SCM_UB_vector(i,j);
            SCM_UB_vectors_out << SCM_UB_vector_ij;
          }
      SCM_UB_vectors_out.close();
    }
}


void RBSCMEvaluation::legacy_read_offline_data_from_files(const std::string & directory_name,
                                                          const bool read_binary_data)
{
  LOG_SCOPE("legacy_read_offline_data_from_files()", "RBSCMEvaluation");

  // The reading mode: DECODE for binary, READ for ASCII
  XdrMODE mode = read_binary_data ? DECODE : READ;

  // The suffix to use for all the files that are written out
  const std::string suffix = read_binary_data ? ".xdr" : ".dat";

  // The string stream we'll use to make the file names
  std::ostringstream file_name;

  // Read in the parameter ranges
  file_name.str("");
  file_name << directory_name << "/parameter_ranges" << suffix;
  std::string continuous_param_file_name = file_name.str();

  // Read in the discrete parameter values
  file_name.str("");
  file_name << directory_name << "/discrete_parameter_values" << suffix;
  std::string discrete_param_file_name = file_name.str();
  read_parameter_data_from_files(continuous_param_file_name,
                                 discrete_param_file_name,
                                 read_binary_data);

  // Read in the bounding box min values
  // Note that there are Q_a values
  file_name.str("");
  file_name << directory_name << "/B_min" << suffix;
  Xdr B_min_in(file_name.str(), mode);

  B_min.clear();
  for (unsigned int i=0; i<rb_theta_expansion->get_n_A_terms(); i++)
    {
      Real B_min_val;
      B_min_in >> B_min_val;
      B_min.push_back(B_min_val);
    }
  B_min_in.close();


  // Read in the bounding box max values
  // Note that there are Q_a values
  file_name.str("");
  file_name << directory_name << "/B_max" << suffix;
  Xdr B_max_in(file_name.str(), mode);

  B_max.clear();
  for (unsigned int i=0; i<rb_theta_expansion->get_n_A_terms(); i++)
    {
      Real B_max_val;
      B_max_in >> B_max_val;
      B_max.push_back(B_max_val);
    }

  // Read in the length of the C_J data
  file_name.str("");
  file_name << directory_name << "/C_J_length" << suffix;
  Xdr C_J_length_in(file_name.str(), mode);

  unsigned int C_J_length;
  C_J_length_in >> C_J_length;
  C_J_length_in.close();

  // Read in C_J_stability_vector
  file_name.str("");
  file_name << directory_name << "/C_J_stability_vector" << suffix;
  Xdr C_J_stability_vector_in(file_name.str(), mode);

  C_J_stability_vector.clear();
  for (unsigned int i=0; i<C_J_length; i++)
    {
      Real C_J_stability_val;
      C_J_stability_vector_in >> C_J_stability_val;
      C_J_stability_vector.push_back(C_J_stability_val);
    }
  C_J_stability_vector_in.close();

  // Read in C_J
  file_name.str("");
  file_name << directory_name << "/C_J" << suffix;
  Xdr C_J_in(file_name.str(), mode);

  // Resize C_J based on C_J_stability_vector and Q_a
  C_J.resize( C_J_length );
  for (auto & params : C_J)
    for (const auto & pr : get_parameters())
      {
        const std::string & param_name = pr.first;
        Real param_value;
        C_J_in >> param_value;
        params.set_value(param_name, param_value);
      }
  C_J_in.close();


  // Read in SCM_UB_vectors get_SCM_UB_vector
  file_name.str("");
  file_name << directory_name << "/SCM_UB_vectors" << suffix;
  Xdr SCM_UB_vectors_in(file_name.str(), mode);

  // Resize SCM_UB_vectors based on C_J_stability_vector and Q_a
  SCM_UB_vectors.resize( C_J_stability_vector.size() );
  for (auto i : index_range(SCM_UB_vectors))
    {
      SCM_UB_vectors[i].resize( rb_theta_expansion->get_n_A_terms() );
      for (unsigned int j=0; j<rb_theta_expansion->get_n_A_terms(); j++)
        {
          SCM_UB_vectors_in >> SCM_UB_vectors[i][j];
        }
    }
  SCM_UB_vectors_in.close();
}

} // namespace libMesh

  _subset_solve_mode(SUBSET_ZERO)
{
  if (this->n_processors() == 1)
    this->_preconditioner_type = ILU_PRECOND;
  else
    this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
}

      if (!this->_preconditioner)
        {
          if (this->n_processors() == 1)
            this->_preconditioner_type = ILU_PRECOND;
          else
            this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
        }

 * restricted here into PETSc's reduced \p r. Otherwise \p R wraps \p r, and the constraint rows
 * are enforced in it
 */
void libmesh_petsc_snes_finish_residual(ResidualContext & rc,
                                        Vec x,
                                        Vec r,
                                        NumericVector<Number> & R)
{
  R.close();

  if (rc.solver->eliminating_constraints())
    {
      rc.solver->restrict_vector(R, r);
      return;
    }

  // Synchronize PETSc x to local solution since the local solution may be changed due to the constraints
  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(rc.sys.solution.get());
  PetscVector<Number> X_global(x, rc.sys.comm());

  X_global.swap(X_sys);
  rc.sys.update();
  X_global.swap(X_sys);

  if (rc.solver->_exact_constraint_enforcement)
    {
      rc.sys.get_dof_map().enforce_constraints_on_residual(rc.sys, &R, rc.sys.current_local_solution.get(), /*homogeneous=*/false);
      R.close();
    }
}

ResidualContext
libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx)

  NonlinearImplicitSystem & sys = solver->system();

  solver->localize_solution(x);

  return ResidualContext(solver, sys);
}

    libmesh_parallel_only(rc.sys.comm());

    libmesh_assert(r);
    PetscVector<Number> R_petsc(r, rc.sys.comm());
    NumericVector<Number> & R =
      rc.solver->eliminating_constraints() ? *rc.sys.rhs : static_cast<NumericVector<Number> &>(R_petsc);

    if (rc.solver->_zero_out_residual)
      R.zero();

          *rc.sys.current_local_solution.get(), &R, &jac, rc.sys);

      jac.close();
      if (rc.solver->_exact_constraint_enforcement && !rc.solver->eliminating_constraints())
        {
          rc.sys.get_dof_map().enforce_constraints_on_jacobian(rc.sys, &jac);
          jac.close();

      libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");


    libmesh_petsc_snes_finish_residual(rc, x, r, R);

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

    libmesh_parallel_only(rc.sys.comm());

    libmesh_assert(r);
    PetscVector<Number> R_petsc(r, rc.sys.comm());
    NumericVector<Number> & R =
      rc.solver->eliminating_constraints() ? *rc.sys.rhs : static_cast<NumericVector<Number> &>(R_petsc);

    if (rc.solver->_zero_out_residual)
      R.zero();

    else
      libmesh_error_msg("Error! Unable to compute residual for forming finite difference Jacobian!");

    libmesh_petsc_snes_finish_residual(rc, x, r, R);

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

    libmesh_parallel_only(rc.sys.comm());

    libmesh_assert(r);
    PetscVector<Number> R_petsc(r, rc.sys.comm());
    NumericVector<Number> & R =
      rc.solver->eliminating_constraints() ? *rc.sys.rhs : static_cast<NumericVector<Number> &>(R_petsc);

    if (rc.solver->_zero_out_residual)
      R.zero();

      libmesh_error_msg("Error! Unable to compute residual for forming finite differenced"
                        "Jacobian-vector products!");

    libmesh_petsc_snes_finish_residual(rc, x, r, R);

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

        Jac->close();

      if (pc && (p_is_shell == PETSC_TRUE))
        PC->close();

      PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
    }

    PC->attach_dof_map(sys.get_dof_map());
    Jac->attach_dof_map(sys.get_dof_map());

    solver->localize_solution(x);

    if (solver->_zero_out_jacobian)
      PC->zero();

    PC->close();
    // A reduced Jacobian has no constraint rows to enforce
    if (solver->_exact_constraint_enforcement && !solver->eliminating_constraints())
      {
        sys.get_dof_map().enforce_constraints_on_jacobian(sys, PC);
        PC->close();

    if (Jac != PC)
      {
        // Assume that shells know what they're doing
        libmesh_assert(!solver->_exact_constraint_enforcement || solver->eliminating_constraints() ||
                       (j_is_mffd == PETSC_TRUE) || (j_is_shell == PETSC_TRUE));
        Jac->close();
      }

      changed_search_direction = false,
      changed_new_soln = false;

    auto call_postcheck = [&](const NumericVector<Number> & old_soln,
                              NumericVector<Number> & search_direction,
                              NumericVector<Number> & new_soln) {
      if (solver->postcheck)
        solver->postcheck(old_soln,
                          search_direction,
                          new_soln,
                          changed_search_direction,
                          changed_new_soln,
                          sys);
      else
        solver->postcheck_object->postcheck(old_soln,
                                            search_direction,
                                            new_soln,
                                            changed_search_direction,
                                            changed_new_soln,
                                            sys);
    };

    if (solver->eliminating_constraints())
      {
        // The user sees full iterates C x + h and C w + h and the full step C y. Restricting the
        // user's changes back to the reduced space keeps each change to an unconstrained entry
        // and projects the iterate back onto the constraints
        auto old_soln = sys.solution->zero_clone();
        auto search_direction = sys.solution->zero_clone();
        auto new_soln = sys.solution->zero_clone();
        solver->prolong(x, *old_soln, /*homogeneous=*/false);
        solver->prolong(y, *search_direction, /*homogeneous=*/true);
        solver->prolong(w, *new_soln, /*homogeneous=*/false);

        call_postcheck(*old_soln, *search_direction, *new_soln);

        if (changed_search_direction)
          solver->restrict_vector(*search_direction, y);
        if (changed_new_soln)
          solver->restrict_vector(*new_soln, w);
      }
    else
      {
        PetscVector<Number> petsc_x(x, sys.comm());
        PetscVector<Number> petsc_y(y, sys.comm());
        PetscVector<Number> petsc_w(w, sys.comm());
        call_postcheck(petsc_x, petsc_y, petsc_w);
      }

    // Record whether the user changed the solution or the search direction.

    auto & sys = solver->system();

    // The precheck solution is the system's current_local_solution, which holds the full iterate
    solver->localize_solution(X);
    libmesh_assert(sys.current_local_solution.get());
    auto & local_soln = *sys.current_local_solution.get();

    if (solver->eliminating_constraints())
      {
        // The user sees the full search direction C Y, and a change to it is restricted back
        auto search_direction = sys.solution->zero_clone();
        solver->prolong(Y, *search_direction, /*homogeneous=*/true);
        solver->precheck_object->precheck(local_soln, *search_direction, petsc_changed, sys);
        if (petsc_changed)
          solver->restrict_vector(*search_direction, Y);
      }
    else
      {
        PetscVector<Number> petsc_y(Y, sys.comm());
        solver->precheck_object->precheck(local_soln, petsc_y, petsc_changed, sys);
      }

    // Record whether the user changed the solution or the search direction.
    if (petsc_changed)

      if (name)
        {
          libmesh_assert(std::string(name).front() != '-');
          libmesh_assert(std::string(name).back() == '_');
          LibmeshPetscCall(SNESSetOptionsPrefix(_snes, name));
        }

      // Attaching a DM to SNES.

        LibmeshPetscCall(DMlibMeshSetSystem(dm, this->system()));

        if (name)
          LibmeshPetscCall(DMSetOptionsPrefix(dm, name));

        LibmeshPetscCall(DMSetFromOptions(dm));
        LibmeshPetscCall(DMSetUp(dm));

template <typename T>
void
PetscNonlinearSolver<T>::build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
                                              void (*computeSubspace)(std::vector<NumericVector<Number> *> &, sys_type &),
                                              MatNullSpace * msp)
{
  parallel_object_only();

  std::vector<NumericVector<Number> *> sp;
  if (computeSubspaceObject)
    (*computeSubspaceObject)(sp, this->system());
  else
    (*computeSubspace)(sp, this->system());

  *msp = LIBMESH_PETSC_NULLPTR;
  if (sp.size())
    {
      PetscInt nmodes = cast_int<PetscInt>(sp.size());

      std::vector<Vec> modes(nmodes);
      std::vector<PetscScalar> dots(nmodes);

      for (PetscInt i=0; i<nmodes; ++i)
        {
          auto pv = cast_ptr<PetscVector<T> *>(sp[i]);

          if (this->eliminating_constraints())
            {
              // A full mode v = C y of the constrained operator gives the reduced mode y, the
              // unconstrained entries of v
              const auto & reduced_map = this->system().get_dof_map().get_constraint_elimination();
              LibmeshPetscCall(VecCreate(this->comm().get(), &modes[i]));
              LibmeshPetscCall(VecSetSizes(modes[i],
                                           cast_int<PetscInt>(reduced_map.n_local_dofs()),
                                           cast_int<PetscInt>(reduced_map.n_dofs())));
              VecType vec_type;
              LibmeshPetscCall(VecGetType(pv->vec(), &vec_type));
              LibmeshPetscCall(VecSetType(modes[i], vec_type));
              LibmeshPetscCall(VecSetUp(modes[i]));
              this->restrict_vector(*pv, modes[i]);
            }
          else
            {
              LibmeshPetscCall(VecDuplicate(pv->vec(), &modes[i]));

              LibmeshPetscCall(VecCopy(pv->vec(), modes[i]));
            }
        }

      // Normalize.
      LibmeshPetscCall(VecNormalize(modes[0], LIBMESH_PETSC_NULLPTR));

      for (PetscInt i=1; i<nmodes; i++)
        {
          // Orthonormalize vec[i] against vec[0:i-1]
          LibmeshPetscCall(VecMDot(modes[i], i, modes.data(), dots.data()));

          LibmeshPetscCall(VecNormalize(modes[i], LIBMESH_PETSC_NULLPTR));
        }

      LibmeshPetscCall(MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes.data(), msp));

      for (PetscInt i=0; i<nmodes; ++i)
        LibmeshPetscCall(VecDestroy(&modes[i]));
    }
}

template <typename T>
std::pair<unsigned int, Real>

  // When eliminating constraints, SNES iterates on the reduced vector y, from which the full
  // solution x = C y + h is formed for every evaluation, and its residual is the reduced one.
  // The full x_in and r_in remain the system's assembly vectors
  WrappedPetsc<Vec> reduced_x, reduced_r;
  this->setup_constraint_elimination();
  if (this->eliminating_constraints())
    {
      const auto & reduced_map = this->system().get_dof_map().get_constraint_elimination();
      const auto n = cast_int<PetscInt>(reduced_map.n_dofs());
      const auto n_local = cast_int<PetscInt>(reduced_map.n_local_dofs());
      LibmeshPetscCall(VecCreate(this->comm().get(), reduced_x.get()));
      LibmeshPetscCall(VecSetSizes(reduced_x, n_local, n));
      VecType vec_type;
      LibmeshPetscCall(VecGetType(x->vec(), &vec_type));
      LibmeshPetscCall(VecSetType(reduced_x, vec_type));
      LibmeshPetscCall(VecSetUp(reduced_x));
      LibmeshPetscCall(VecDuplicate(reduced_x, reduced_r.get()));

      // The initial guess is the unconstrained part of the current solution
      x->close();
      this->restrict_vector(*x, reduced_x);
    }
  const Vec snes_x = this->eliminating_constraints() ? static_cast<Vec>(reduced_x) : x->vec();
  const Vec snes_r = this->eliminating_constraints() ? static_cast<Vec>(reduced_r) : r->vec();

  PetscInt n_iterations =0;
  // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal

        }
    }

  LibmeshPetscCall(SNESSetFunction (_snes, snes_r, libmesh_petsc_snes_residual, this));

  // Only set the jacobian function if we've been provided with something to call.
  // This allows a user to set their own jacobian function if they want to

  // (petsc/petsc@154060b), so this code block should be safe to use
  // in 3.5.0 and later.
#if !PETSC_VERSION_LESS_THAN(3,6,0)
  LibmeshPetscCall(SNESSetSolution(_snes, snes_x));
#endif
  LibmeshPetscCall(SNESSetUp(_snes));

  // Only set the nullspace if we have a way of computing it and the result is non-empty.
  if (this->nullspace || this->nullspace_object)
    {
      WrappedPetsc<MatNullSpace> msp;
      this->build_mat_null_space(this->nullspace_object, this->nullspace, msp.get());
      if (msp)
        {
          LibmeshPetscCall(MatSetNullSpace(J, msp));
          if (P != J)
            LibmeshPetscCall(MatSetNullSpace(P, msp));
        }
    }

#endif
  }

  LibmeshPetscCall(SNESSolve (_snes, LIBMESH_PETSC_NULLPTR, snes_x));

  if (this->eliminating_constraints())
    this->prolong(snes_x, *x, /*homogeneous=*/false);

  LibmeshPetscCall(SNESGetIterationNumber(_snes, &n_iterations));

}

template <typename T>
void PetscNonlinearSolver<T>::setup_constraint_elimination()
{
  _retained_is.destroy();

  const DofMap & dof_map = this->system().get_dof_map();
  if (!dof_map.has_constraint_elimination())
    return;

  const auto & retained = dof_map.get_constraint_elimination().local_retained_dofs();
  std::vector<PetscInt> indices(retained.begin(), retained.end());
  LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
                                   cast_int<PetscInt>(indices.size()),
                                   indices.data(),
                                   PETSC_COPY_VALUES,

}

template <typename T>
void PetscNonlinearSolver<T>::localize_solution(Vec x)
{
  sys_type & sys = this->system();

  if (this->eliminating_constraints())
    {
      // x is the reduced iterate y; the system solution becomes C y + h
      this->prolong(x, *sys.solution, /*homogeneous=*/false);
      sys.update();
      return;
    }

  PetscVector<T> & X_sys = *cast_ptr<PetscVector<T> *>(sys.solution.get());
  PetscVector<T> X_global(x, sys.comm());

  // Use the system's update() to get a good local version of the
  // parallel solution.  This operation does not modify the incoming
  // "x" vector, it only localizes information from "x" into
  // sys.current_local_solution.
  X_global.swap(X_sys);
  sys.update();
  X_global.swap(X_sys);

  // Enforce constraints (if any) exactly on the
  // current_local_solution.  This is the solution vector that is
  // actually used in residual and Jacobian evaluations, and is
  // not locked by debug-enabled PETSc the way that "x" is.
  if (this->_exact_constraint_enforcement)
    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
}

template <typename T>
void PetscNonlinearSolver<T>::prolong(Vec y, NumericVector<T> & x, const bool homogeneous) const
{
  libmesh_assert(this->eliminating_constraints());
  auto & petsc_x = cast_ref<PetscVector<T> &>(x);
  libmesh_assert_equal_to(petsc_x.type(), PARALLEL);

  // The constrained entries of x are overwritten below, so their prior values do not matter
  LibmeshPetscCall(VecISCopy(petsc_x.vec(), _retained_is, SCATTER_FORWARD, y));
  x.close();
  this->system().get_dof_map().enforce_constraints_exactly(this->system(), &x, homogeneous);
}

template <typename T>
void PetscNonlinearSolver<T>::restrict_vector(const NumericVector<T> & x, Vec y) const
{
  libmesh_assert(this->eliminating_constraints());
  const auto & petsc_x = cast_ref<const PetscVector<T> &>(x);
  LibmeshPetscCall(VecISCopy(petsc_x.vec(), _retained_is, SCATTER_REVERSE, y));
}

//------------------------------------------------------------------
// Explicit instantiations

  std::map<unsigned int, std::string> * dnames = LIBMESH_PETSC_NULLPTR,
                                      * enames = LIBMESH_PETSC_NULLPTR;
  if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
    name += ":dec:var:";
    dnames = dlm->varnames;
  }
  if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
    name += ":dec:block:";
    dnames = dlm->blocknames;
  }
  if (dnames) {
    for (auto decomp : *dlm->decomposition) {
      for (std::set<unsigned int>::iterator dit_begin = decomp.begin(),
                                            dit = dit_begin,
                                            dit_end   = decomp.end();
           dit != dit_end; ++dit) {
        unsigned int id = *dit;
        if (dit != dit_begin)
          name += ",";
        name += (*dnames)[id];
      }
      name += ";";
    }
  }
  if (dlm->embedding_type == DMLIBMESH_FIELD_EMBEDDING) {

 * variable \p var on \p elem: reduced indices when \p sys eliminates constrained degrees of
 * freedom, and full ones otherwise
 */
static void DMlibMeshInsertLocalDofs_Private(const NonlinearImplicitSystem & sys,
                                             const Elem * elem,
                                             unsigned int var,
                                             std::set<numeric_index_type> & indices)
{
  const DofMap & dof_map = sys.get_dof_map();
  std::vector<dof_id_type> evindices;
  dof_map.dof_indices(elem, evindices, var);

  const ConstraintEliminationDofMap * reduced =
    dof_map.has_constraint_elimination() ? &dof_map.get_constraint_elimination() : nullptr;

  for (const auto dof : evindices)
    if (dof_map.local_index(dof))
      {
        if (!reduced)
          indices.insert(dof);
        else if (const auto r = reduced->find_reduced_dof(dof); r != DofObject::invalid_id)
          indices.insert(r);
      }
}

#undef __FUNCT__
#define __FUNCT__ "DMCreateFieldDecomposition_libMesh"
static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
{
  PetscFunctionBegin;
  DM_libMesh     * dlm = (DM_libMesh *)(dm->data);
  NonlinearImplicitSystem * sys = dlm->sys;
  IS emb;
  if (dlm->decomposition_type != DMLIBMESH_FIELD_DECOMPOSITION) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);

  *len = cast_int<unsigned int>(dlm->decomposition->size());
  if (namelist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(char *), namelist));}
  if (islist)   {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(IS),    islist));}
  if (dmlist)   {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(DM),    dmlist));}
  for (auto d : index_range(*dlm->decomposition)) {
    std::set<numeric_index_type>         dindices;
    std::string                          dname;
    std::map<std::string, unsigned int>  dvarids;
    std::map<unsigned int, std::string>  dvarnames;
    unsigned int                         dvcount = 0;
    for (const auto & v : (*dlm->decomposition)[d]) {
      std::string vname = (*dlm->varnames)[v];
      dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
      dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
      if (!dvcount) dname = vname;
      else   dname += "_" + vname;
      ++dvcount;
      if (!islist) continue;
      // Iterate only over this DM's blocks.
      for (const auto & pr : *(dlm->blockids)) {
        const subdomain_id_type sbd_id = cast_int<subdomain_id_type>(pr.second);
        for (const auto & elem :
               as_range(sys->get_mesh().active_local_subdomain_elements_begin(sbd_id),
                        sys->get_mesh().active_local_subdomain_elements_end(sbd_id))) {
          DMlibMeshInsertLocalDofs_Private(*sys, elem, v, dindices);
        }
      }
    }
    if (namelist) {
      LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(),(*namelist)+d));
    }
    if (islist) {
      IS dis;
      PetscInt * darray;
      LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray));
      numeric_index_type i = 0;
      for (const auto & id : dindices) {
        darray[i] = id;
        ++i;
      }
      LibmeshPetscCallQ(ISCreateGeneral(((PetscObject)dm)->comm,
                                        cast_int<PetscInt>(dindices.size()),
                                        darray, PETSC_OWN_POINTER, &dis));
      if (dlm->embedding) {
        /* Create a relative embedding into the parent's index space. */
        LibmeshPetscCallQ(ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb));
        PetscInt elen, dlen;

        dis = emb;
      }
      else {
        emb = dis;
      }
      (*islist)[d] = dis;
    }
    if (dmlist) {
      DM ddm;
      LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, &ddm));
      LibmeshPetscCallQ(DMSetType(ddm, DMLIBMESH));

      (*dmlist)[d] = ddm;
    }
  }
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}

#undef __FUNCT__

      for ( ; el != end_el; ++el) {
        const Elem * elem = *el;
        // Iterate only over this DM's variables.
        for (const auto & pr : *(dlm->varids))
          DMlibMeshInsertLocalDofs_Private(*sys, elem, pr.second, dindices);
      }
    }
    if (namelist) {

#undef  __FUNCT__
#define __FUNCT__ "DMlibMeshCreateFieldDecompositionDM"
PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dvarlists, DM * ddm)
{
  PetscBool islibmesh;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
  if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
  if (dnumber < 0) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %" LIBMESH_PETSCINT_FMT " of decomposition parts", dnumber);
#if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
  PetscAssertPointer(ddm,5);
#else
  PetscValidPointer(ddm,5);
#endif
  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
  LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, ddm));
  LibmeshPetscCallQ(DMSetType(*ddm, DMLIBMESH));
  DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
  ddlm->sys = dlm->sys;
  // Each DM owns, and DMDestroy_libMesh() deletes, its own variable and block maps
  *ddlm->varids = *dlm->varids;
  *ddlm->varnames = *dlm->varnames;
  *ddlm->blockids = *dlm->blockids;
  *ddlm->blocknames = *dlm->blocknames;
  ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
  ddlm->decomposition_type = DMLIBMESH_FIELD_DECOMPOSITION;
  if (dnumber) {
    for (PetscInt d = 0; d < dnumber; ++d) {
      if (dsizes[d] < 0) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %" LIBMESH_PETSCINT_FMT " of decomposition part %" LIBMESH_PETSCINT_FMT, dsizes[d],d);
      ddlm->decomposition->push_back(std::set<unsigned int>());

    }
  }
  else { // Empty splits indicate default: split all variables with one per split.
    PetscInt d = 0;
    for (const auto & pr : (*ddlm->varids)) {
      ddlm->decomposition->push_back(std::set<unsigned int>());
      unsigned int vid = pr.second;
      (*ddlm->decomposition)[d].insert(vid);
      ++d;
    }
  }
  LibmeshPetscCallQ(DMlibMeshSetUpName_Private(*ddm));
  LibmeshPetscCallQ(DMSetFromOptions(*ddm));
  LibmeshPetscCallQ(DMSetUp(*ddm));
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}

#undef  __FUNCT__

  DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
  ddlm->sys = dlm->sys;
  // Each DM owns, and DMDestroy_libMesh() deletes, its own variable and block maps
  *ddlm->varids   = *dlm->varids;
  *ddlm->varnames = *dlm->varnames;
  *ddlm->blockids   = *dlm->blockids;
  *ddlm->blocknames = *dlm->blocknames;
  ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
  ddlm->decomposition_type = DMLIBMESH_DOMAIN_DECOMPOSITION;
  if (dnumber) {

  NonlinearImplicitSystem * _sys;
  LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
  NonlinearImplicitSystem & sys = *_sys;
  if (sys.has_constraint_elimination())
    SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_SUP, "DMlibMesh does not evaluate residuals of systems that eliminate constraints; PetscNonlinearSolver does");
  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
  PetscVector<Number> & R_sys = *cast_ptr<PetscVector<Number> *>(sys.rhs);
  PetscVector<Number> X_global(x, _sys->comm()), R(r, _sys->comm());

  NonlinearImplicitSystem * _sys;
  LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
  NonlinearImplicitSystem & sys = *_sys;
  if (sys.has_constraint_elimination())
    SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_SUP, "DMlibMesh does not evaluate Jacobians of systems that eliminate constraints; PetscNonlinearSolver does");

  libmesh_assert(pc);
  libmesh_assert(jac);

  // Users compute bounds on full vectors. When eliminating constraints, SNES's xl and xu are
  // reduced, and they take the bounds of the unconstrained degrees of freedom
  const bool reduced = sys.has_constraint_elimination();
  std::unique_ptr<NumericVector<Number>> full_l, full_u;
  if (reduced)
    {
      full_l = sys.solution->zero_clone();
      full_u = sys.solution->zero_clone();
    }
  PetscVector<Number> XL(xl, sys.comm()), XU(xu, sys.comm());
  NumericVector<Number> & L = reduced ? *full_l : static_cast<NumericVector<Number> &>(XL);
  NumericVector<Number> & U = reduced ? *full_u : static_cast<NumericVector<Number> &>(XU);
  L.add(-petsc_inf);
  U.add(petsc_inf);
  L.close();
  U.close();

  if (sys.nonlinear_solver->bounds != nullptr)
    sys.nonlinear_solver->bounds(L,U,sys);
  else if (sys.nonlinear_solver->bounds_object != nullptr)
    sys.nonlinear_solver->bounds_object->bounds(L,U, sys);
  else
    SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this libMesh object");

  if (reduced)
    {
      auto & solver = cast_ref<PetscNonlinearSolver<Number> &>(*sys.nonlinear_solver);
      L.close();
      U.close();
      solver.restrict_vector(L, xl);
      solver.restrict_vector(U, xu);
    }

  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);

  NumericVector<Number> * nv = (dlm->sys->solution).get();
  PetscVector<Number> *   pv = dynamic_cast<PetscVector<Number> *>(nv);
  Vec                    v  = pv->vec();
  const DofMap & dof_map = dlm->sys->get_dof_map();
  /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves aren't going to be easily available.
     Should work fine for getting vectors out for linear subproblem solvers. */
  if (dlm->embedding || dof_map.has_constraint_elimination()) {
    /* The size of an embedded subproblem, or of the reduced space SNES solves on */
    PetscInt n;
    LibmeshPetscCallQ(VecCreate(((PetscObject)v)->comm, x));
    if (dlm->embedding)
      LibmeshPetscCallQ(ISGetLocalSize(dlm->embedding, &n));
    else
      n = cast_int<PetscInt>(dof_map.get_constraint_elimination().n_local_dofs());
    LibmeshPetscCallQ(VecSetSizes(*x,n,PETSC_DETERMINE));
    LibmeshPetscCallQ(VecSetType(*x,((PetscObject)v)->type_name));
    LibmeshPetscCallQ(VecSetFromOptions(*x));
    LibmeshPetscCallQ(VecSetUp(*x));
  }
  else {
    LibmeshPetscCallQ(VecDuplicate(v,x));

      LibmeshPetscCall(EPSSetType (_eps, EPSARNOLDI)); return;
    case LANCZOS:
      LibmeshPetscCall(EPSSetType (_eps, EPSLANCZOS)); return;
    case KRYLOVSCHUR:
      LibmeshPetscCall(EPSSetType (_eps, EPSKRYLOVSCHUR)); return;
      // case ARPACK:
      // LibmeshPetscCall(EPSSetType (_eps, (char *) EPSARPACK)); return;

        LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_MAGNITUDE));
        return;
      }
    case LARGEST_REAL:
      {
        LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_LARGEST_REAL));
        return;
      }
    case SMALLEST_REAL:
      {
        LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_REAL));
        return;
      }
    case LARGEST_IMAGINARY:

  //---------------------------------------------------------------
  // This function is called by Tao to evaluate the equality constraints at x
  PetscErrorCode
  __libmesh_tao_equality_constraints(Tao /*tao*/, Vec x, Vec ce, void * ctx)
  {
    PetscFunctionBegin;

    // We'll use current_local_solution below, so let's ensure that it's consistent
    // with the vector x that was passed in.
    PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
    PetscVector<Number> X(x, sys.comm());

    // Perform a swap so that sys.solution points to the input vector
    // "x", update sys.current_local_solution based on "x", then swap
    // back.
    X.swap(X_sys);
    sys.update();
    X.swap(X_sys);

    // We'll also pass the constraints vector ce into the assembly routine
    // so let's make a PETSc vector for that too.
    PetscVector<Number> eq_constraints(ce, sys.comm());

    // Clear the gradient prior to assembly
    eq_constraints.zero();

    // Enforce constraints exactly on the current_local_solution.
    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());

    if (solver->equality_constraints_object != nullptr)
      solver->equality_constraints_object->equality_constraints(*(sys.current_local_solution), eq_constraints, sys);
    else
      libmesh_error_msg("Constraints function not defined in __libmesh_tao_equality_constraints");

    eq_constraints.close();

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

  //---------------------------------------------------------------
  // This function is called by Tao to evaluate the Jacobian of the
  // equality constraints at x
  PetscErrorCode
  __libmesh_tao_equality_constraints_jacobian(Tao /*tao*/, Vec x, Mat J, Mat Jpre, void * ctx)
  {
    PetscFunctionBegin;

    // We'll use current_local_solution below, so let's ensure that it's consistent
    // with the vector x that was passed in.
    PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
    PetscVector<Number> X(x, sys.comm());

    // Perform a swap so that sys.solution points to the input vector
    // "x", update sys.current_local_solution based on "x", then swap
    // back.
    X.swap(X_sys);
    sys.update();
    X.swap(X_sys);

    // Let's also wrap J and Jpre in PetscMatrix objects for convenience
    PetscMatrix<Number> J_petsc(J, sys.comm());
    PetscMatrix<Number> Jpre_petsc(Jpre, sys.comm());

    // Enforce constraints exactly on the current_local_solution.
    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());

    if (solver->equality_constraints_jacobian_object != nullptr)
      solver->equality_constraints_jacobian_object->equality_constraints_jacobian(*(sys.current_local_solution), J_petsc, sys);
    else
      libmesh_error_msg("Constraints function not defined in __libmesh_tao_equality_constraints_jacobian");

    J_petsc.close();
    Jpre_petsc.close();

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

  //---------------------------------------------------------------
  // This function is called by Tao to evaluate the inequality constraints at x
  PetscErrorCode
  __libmesh_tao_inequality_constraints(Tao /*tao*/, Vec x, Vec cineq, void * ctx)
  {
    PetscFunctionBegin;

    // We'll use current_local_solution below, so let's ensure that it's consistent
    // with the vector x that was passed in.
    PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
    PetscVector<Number> X(x, sys.comm());

    // Perform a swap so that sys.solution points to the input vector
    // "x", update sys.current_local_solution based on "x", then swap
    // back.
    X.swap(X_sys);
    sys.update();
    X.swap(X_sys);

    // We'll also pass the constraints vector ce into the assembly routine
    // so let's make a PETSc vector for that too.
    PetscVector<Number> ineq_constraints(cineq, sys.comm());

    // Clear the gradient prior to assembly
    ineq_constraints.zero();

    // Enforce constraints exactly on the current_local_solution.
    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());

    if (solver->inequality_constraints_object != nullptr)
      solver->inequality_constraints_object->inequality_constraints(*(sys.current_local_solution), ineq_constraints, sys);
    else
      libmesh_error_msg("Constraints function not defined in __libmesh_tao_inequality_constraints");

    ineq_constraints.close();

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

  //---------------------------------------------------------------
  // This function is called by Tao to evaluate the Jacobian of the
  // equality constraints at x
  PetscErrorCode
  __libmesh_tao_inequality_constraints_jacobian(Tao /*tao*/, Vec x, Mat J, Mat Jpre, void * ctx)
  {
    PetscFunctionBegin;

    // We'll use current_local_solution below, so let's ensure that it's consistent
    // with the vector x that was passed in.
    PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
    PetscVector<Number> X(x, sys.comm());

    // Perform a swap so that sys.solution points to the input vector
    // "x", update sys.current_local_solution based on "x", then swap
    // back.
    X.swap(X_sys);
    sys.update();
    X.swap(X_sys);

    // Let's also wrap J and Jpre in PetscMatrix objects for convenience
    PetscMatrix<Number> J_petsc(J, sys.comm());
    PetscMatrix<Number> Jpre_petsc(Jpre, sys.comm());

    // Enforce constraints exactly on the current_local_solution.
    sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());

    if (solver->inequality_constraints_jacobian_object != nullptr)
      solver->inequality_constraints_jacobian_object->inequality_constraints_jacobian(*(sys.current_local_solution), J_petsc, sys);
    else
      libmesh_error_msg("Constraints function not defined in __libmesh_tao_inequality_constraints_jacobian");

    J_petsc.close();
    Jpre_petsc.close();

    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }

} // end extern "C"
//---------------------------------------------------------------------

  if (this->lower_and_upper_bounds_object)
    {
      // Need to actually compute the bounds vectors first
      this->lower_and_upper_bounds_object->lower_and_upper_bounds(this->system());

      LibmeshPetscCall(TaoSetVariableBounds(_tao,
                                            lb->vec(),
                                            ub->vec()));
    }

  if (this->equality_constraints_object)
    LibmeshPetscCall(TaoSetEqualityConstraintsRoutine(_tao, ceq->vec(), __libmesh_tao_equality_constraints, this));

  if (this->equality_constraints_jacobian_object)
    LibmeshPetscCall(TaoSetJacobianEqualityRoutine(_tao,
                                                   ceq_jac->mat(),
                                                   ceq_jac->mat(),
                                                   __libmesh_tao_equality_constraints_jacobian,

  // Optionally set inequality constraints
  if (this->inequality_constraints_object)
    LibmeshPetscCall(TaoSetInequalityConstraintsRoutine(_tao, cineq->vec(), __libmesh_tao_inequality_constraints, this));

  // Optionally set inequality constraints Jacobian
  if (this->inequality_constraints_jacobian_object)
    LibmeshPetscCall(TaoSetJacobianInequalityRoutine(_tao,
                                                     cineq_jac->mat(),
                                                     cineq_jac->mat(),
                                                     __libmesh_tao_inequality_constraints_jacobian,

template <typename T>
void TaoOptimizationSolver<T>::get_dual_variables()
{
  LOG_SCOPE("get_dual_variables()", "TaoOptimizationSolver");

  PetscVector<T> * lambda_ineq_petsc =
    cast_ptr<PetscVector<T> *>(this->system().lambda_ineq.get());

  Vec lambda_eq_petsc_vec = lambda_eq_petsc->vec();
  Vec lambda_ineq_petsc_vec = lambda_ineq_petsc->vec();

  LibmeshPetscCall(TaoGetDualVariables(_tao,
                                       &lambda_eq_petsc_vec,
                                       &lambda_ineq_petsc_vec));
}


template <typename T>

  // Now erase the condensed dofs
  for (const auto dof : global_dirichlet_dofs_set)
    if ((dof_map.first_dof() <= dof) && (dof < dof_map.end_dof()))
      local_non_condensed_dofs_set.erase(dof);

  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()

  if (this->get_mesh_x_var() != libMesh::invalid_uint)
    for (unsigned int i=0; i != n_nodes; ++i)
      const_cast<Elem &>(this->get_elem()).point(i)(0) =
        libmesh_real(this->get_elem_solution(this->get_mesh_x_var())(i));

  if (this->get_mesh_y_var() != libMesh::invalid_uint)
    for (unsigned int i=0; i != n_nodes; ++i)
      const_cast<Elem &>(this->get_elem()).point(i)(1) =
        libmesh_real(this->get_elem_solution(this->get_mesh_y_var())(i));

  if (this->get_mesh_z_var() != libMesh::invalid_uint)
    for (unsigned int i=0; i != n_nodes; ++i)
      const_cast<Elem &>(this->get_elem()).point(i)(2) =
        libmesh_real(this->get_elem_solution(this->get_mesh_z_var())(i));
  //    }
  // FIXME - If the coordinate data is not in our own system, someone
  // had better get around to implementing that... - RHS

void NonlinearImplicitSystem::create_static_condensation()
{
  Parent::create_static_condensation();
  this->setup_static_condensation_preconditioner(*nonlinear_solver);
}

}


void OptimizationSystem::
initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
{
  numeric_index_type n_eq_constraints =

  // Assign rows to each processor as evenly as possible
  unsigned int n_procs = comm().size();
  numeric_index_type n_local_rows = n_eq_constraints / n_procs;
  if (comm().rank() < (n_eq_constraints % n_procs))
    n_local_rows++;

  C_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
  lambda_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);

  // Get the maximum number of non-zeros per row
  numeric_index_type max_nnz = 0;
  for (numeric_index_type i=0; i<n_eq_constraints; i++)
    {
      numeric_index_type nnz =
        cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());

        max_nnz = nnz;
    }

  C_eq_jac->init(n_eq_constraints,
                 get_dof_map().n_dofs(),
                 n_local_rows,
                 get_dof_map().n_local_dofs(),
                 max_nnz,
                 max_nnz);

  eq_constraint_jac_sparsity = constraint_jac_sparsity;
}


void OptimizationSystem::
initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
{
  numeric_index_type n_ineq_constraints =

  // Assign rows to each processor as evenly as possible
  unsigned int n_procs = comm().size();
  numeric_index_type n_local_rows = n_ineq_constraints / n_procs;
  if (comm().rank() < (n_ineq_constraints % n_procs))
    n_local_rows++;

  C_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
  lambda_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);

  // Get the maximum number of non-zeros per row
  numeric_index_type max_nnz = 0;
  for (numeric_index_type i=0; i<n_ineq_constraints; i++)
    {
      numeric_index_type nnz =
        cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());

        max_nnz = nnz;
    }

  C_ineq_jac->init(n_ineq_constraints,
                   get_dof_map().n_dofs(),
                   n_local_rows,
                   get_dof_map().n_local_dofs(),
                   max_nnz,
                   max_nnz);

  ineq_constraint_jac_sparsity = constraint_jac_sparsity;
}


void OptimizationSystem::solve ()

    this->prefix_with_name(true);
  if (libMesh::on_command_line("--" + name_in + "-static-condensation"))
    this->create_static_condensation();
  if (libMesh::on_command_line("--" + name_in + "-eliminate-constraints"))
    this->create_constraint_elimination();
}

  // But our "basic"/"null" state may still have a StaticCondensation
  if (libMesh::on_command_line("--" + _sys_name + "-static-condensation"))
    this->create_static_condensation();
  if (libMesh::on_command_line("--" + _sys_name + "-eliminate-constraints"))
    this->create_constraint_elimination();
}

  // Now finally after dof distribution and construction of any
  // possible constraints, we may init any constraint elimination
  // data and then the static condensation data nested in it
  _dof_map->reinit_constraint_elimination();
  _dof_map->reinit_static_condensation();

  // Resize the solution conformal to the current mesh

  // Make sure our constraint elimination and static condensation dof maps are up-to-date before
  // we init any matrices on their reduced spaces
  this->get_dof_map().reinit_constraint_elimination();
  this->get_dof_map().reinit_static_condensation();

  if (!_matrices.empty() && !_basic_system_only)

                                                    this->get_dof_map().get_static_condensation());
    }
#ifdef LIBMESH_HAVE_PETSC
  else if (this->get_dof_map().has_constraint_elimination())
    {
      if (mat_build_type == MatrixBuildType::DIAGONAL)
        libmesh_error_msg(
            "We do not currently support constraint elimination for the diagonal matrix type");
      matrix = std::make_unique<ReducedSparseMatrix<Number>>(
          this->comm(), this->get_dof_map().get_constraint_elimination());
    }
#endif
  else

  this->get_dof_map().create_static_condensation(this->get_mesh(), *this);
}

void System::create_constraint_elimination()
{
  this->get_dof_map().create_constraint_elimination(this->get_mesh(), *this);
}

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

bool System::has_static_condensation() const

        }
    }

  return nullptr;
}