// (*) 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)

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);
}

   * 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

   * 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);

        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 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);

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

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

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

        }
      else if (elem->point(5) == min_point)
        {
          if (!elem->positive_face_orientation(2))
            {
              // Case 7: 180 degree rotation
              i01 = s0+1+e;

          else
            {
              // Case 8: flip about 1-3 diagonal
              i01 = s1+1+e;
              i2 = s0;
              zeta = 1-2*xe_fraction;
              const Real xe = (1-zeta_saved)/2;
              xi_eta(1) = xe*xe_scale;
              xi_eta(0) = xe_scale - xi_eta(1);
            }
        }
    }

        }
      else if (elem->point(5) == min_point)
        {
          if (!elem->positive_face_orientation(3))
            {
              // Case 3: 2->5->3->0->2 rotation
              i01 = s1+1+2*e;
              i2 = s0;
              zeta = 1-2*xe_fraction;
              const Real xe = (zeta_saved+1)/2;
              xi_eta(1) = xe_scale - xe*xe_scale;
            }
          else
            {

      // 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);
}

  if (_boundary_node_id.empty())
    {
      libMesh::out << "No boundary node IDs have been added: cannot build side list!" << std::endl;
      return;
    }

  // For avoiding extraneous element side construction

}


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));

  // ever.
  if (_elem_type != TRI3 &&
      _elem_type != TRI6 &&
      _elem_type != TRI7)
    libmesh_not_implemented();

  // If we have no explicit segments defined, we may get them from

      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);
}

    case TRI6:
      _mesh.all_second_order();
      break;
    case TRI7:
      _mesh.all_complete_order();
      break;
    default:
      libmesh_not_implemented();
    }

  // Moving boundary mid-edge nodes can displace the TRI7 interior node
  // and tangle the element map, so fix up and verify afterwards.
  if (_elem_type == TRI7)
    this->fixup_tri7_center_nodes();

  this->verify_quadratic_elements();
}


void TriangulatorInterface::fixup_tri7_center_nodes()
{
  libmesh_assert_equal_to(_elem_type, TRI7);

  // Place the interior node at the image of the reference centroid
  // (xi, eta) = (1/3, 1/3) under the curved Tri6 map, using the Tri6

  static const Real wv = -Real(1)/9;
  static const Real wm =  Real(4)/9;

  for (Elem * elem : _mesh.element_ptr_range())
    {
      libmesh_assert_equal_to(elem->n_vertices(), 3);
      libmesh_assert_equal_to(elem->n_nodes(), 7u);

      elem->point(6) = wv * (elem->point(0) +
                             elem->point(1) +
                             elem->point(2)) +
                       wm * (elem->point(3) +
                             elem->point(4) +
                             elem->point(5));
    }
}


void TriangulatorInterface::verify_quadratic_elements()
{
  if (_elem_type != TRI6 && _elem_type != TRI7)
    return;

  // Once fixup_tri7_center_nodes() has placed node 6, the TRI6 and TRI7
  // mappings coincide and this Tri6 formula serves both.

                                      Real(0),   Real(1)/2, Real(1)/2,
                                      Real(1)/3};

  for (Elem * elem : _mesh.element_ptr_range())
    {
      libmesh_assert_equal_to(elem->n_vertices(), 3);
      libmesh_assert_greater_equal(elem->n_nodes(), 6u);

      const Point & x0 = elem->point(0);
      const Point & x1 = elem->point(1);
      const Point & x2 = elem->point(2);
      const Point & x3 = elem->point(3);
      const Point & x4 = elem->point(4);
      const Point & x5 = elem->point(5);

      // Tri6 mapping derivative coefficients (see Tri6::volume()):
      // dx/dxi = xi*a1 + eta*b1 + c1, dx/deta = xi*b1 + eta*b2 + c2.
      const Point a1 =  4*x0 + 4*x1 - 8*x3;
      const Point b1 =  4*x0 - 4*x3 + 4*x4 - 4*x5;
      const Point c1 = -3*x0 - 1*x1 + 4*x3;
      const Point b2 =  4*x0 + 4*x2 - 8*x5;
      const Point c2 = -3*x0 - 1*x2 + 4*x5;

      // Scale the tolerance by the straight-edge triangle area, which
      // is strictly positive for the valid TRI3 poly2tri input.
      const Real ref_area = 0.5 * cross_norm(x1 - x0, x2 - x0);
      const Real jac_tol = TOLERANCE * ref_area;

      Real min_jac = std::numeric_limits<Real>::max();
      unsigned int worst_sample = 0;
      for (unsigned int s = 0; s != 7; ++s)
        {
          const Real xi  = xi_samples[s];
          const Real eta = eta_samples[s];
          const Point dxi  = xi*a1 + eta*b1 + c1;
          const Point deta = xi*b1 + eta*b2 + c2;
          // z-component of the cross product; the elements are planar.
          const Real jac = dxi(0)*deta(1) - dxi(1)*deta(0);
          if (jac < min_jac)
            {
              min_jac = jac;
              worst_sample = s;
            }
        }

      if (min_jac > jac_tol)
        continue;

      // Build a diagnostic naming every snapped boundary side on this
      // element so the user can immediately see which curved-boundary
      // input caused the tangle.
      std::ostringstream sides;
      for (unsigned int n = 0; n != 3; ++n)
        if (!elem->neighbor_ptr(n))
          {
            const Point straight =
              0.5 * (elem->point(n) + elem->point((n+1) % 3));
            sides << " (boundary side " << n
                  << ": straight midpoint " << straight
                  << ", snapped midpoint " << elem->point(n+3) << ")";
          }

      libmesh_error_msg(
        "TriangulatorInterface: snapping a boundary midpoint produced a "
        "tangled quadratic triangle (element " << elem->id()
        << ", non-positive Jacobian " << min_jac

        << " Refine the boundary discretization so that recorded "
        "midpoints lie closer to their straight-line midpoints, "
        "then retry.");
    }
}

          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;

    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)

  //      (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

              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;
        }

      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 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 ()

  // On a parallel mesh we might not yet have a full bounding box
  if (!mesh.is_serial())
    {
      mesh.comm().min(_lower_bound);
      mesh.comm().max(_upper_bound);
    }

  this->fill(mesh);

        }
    }

  return nullptr;
}