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

    // Pack individual integers into vector
    std::vector<int> buffer =
      {num_dim, num_elem, num_elem_blk, num_node_sets,
       num_side_sets, num_elem_sets, num_edge_blk, num_edge,
       num_face_blk, num_face};

    // broadcast integers
    comm.broadcast(buffer);

    num_elem_sets = buffer[ctr++];
    num_edge_blk  = buffer[ctr++];
    num_edge      = buffer[ctr++];
    num_face_blk  = buffer[ctr++];
    num_face      = buffer[ctr++];
  }

  std::vector<char> title;

   * 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
            {

void C0Polyhedron::connectivity(const unsigned int sf,
                                const IOPackage iop,
                                std::vector<dof_id_type> & conn) const
{
  libmesh_assert_less (sf, this->n_sub_elem());
  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);

  const auto & subtet = this->_triangulation[sf];

  switch (iop)
    {
    case TECPLOT:
      {
        conn.resize(8);
        conn[0] = this->node_id(subtet[0])+1;
        conn[1] = this->node_id(subtet[1])+1;
        conn[2] = this->node_id(subtet[2])+1;
        conn[3] = this->node_id(subtet[2])+1;
        conn[4] = this->node_id(subtet[3])+1;
        conn[5] = this->node_id(subtet[3])+1;
        conn[6] = this->node_id(subtet[3])+1;
        conn[7] = this->node_id(subtet[3])+1;
        return;
      }

    case VTK:
      {
        conn.resize(4);
        conn[0] = this->node_id(subtet[0]);
        conn[1] = this->node_id(subtet[1]);
        conn[2] = this->node_id(subtet[2]);
        conn[3] = this->node_id(subtet[3]);
        return;
      }

    default:
      libmesh_error_msg("Unsupported IO package " << iop);
    }
}

}


std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
{
  std::pair<Real, Real> bounds;

  switch (q)
    {
    case EDGE_LENGTH_RATIO:
      bounds.first  = 1.;
      bounds.second = 4.;
      break;

    case MIN_ANGLE:
      bounds.first  = 30.;
      bounds.second = 180.;
      break;

    case MAX_ANGLE:
      bounds.first  = 60.;
      bounds.second = 180.;
      break;

    case JACOBIAN:
    case SCALED_JACOBIAN:
      bounds.first  = 0.5;
      bounds.second = 1.;
      break;

    default:
      libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;

      bounds.second = -1;
    }

  return bounds;
}

      // 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::vector<ElemQuality> Quality::valid(const ElemType t)
{
  std::vector<ElemQuality> v;

  switch (t)
    {
    case EDGE2:
    case EDGE3:

        break;
      }
    case C0POLYGON:
    case C0POLYHEDRON:
      {
        v = {
          EDGE_LENGTH_RATIO,
          JACOBIAN,
          SCALED_JACOBIAN,
          MAX_ANGLE,
          MIN_ANGLE
        };

        break;
      }

      libmesh_error_msg("Undefined element type!");
    }

  return v;
}

} // namespace libMesh

}


std::pair<Real, Real> Polygon::qual_bounds (const ElemQuality q) const
{
  std::pair<Real, Real> bounds;

  switch (q)
    {
    case EDGE_LENGTH_RATIO:
      bounds.first  = 1.;
      bounds.second = 4.;
      break;

    case MIN_ANGLE:
      bounds.first  = 90. - (180./this->n_sides());
      bounds.second = 180. - (360./this->n_sides());
      break;

    case MAX_ANGLE:
      bounds.first  = 180. - (360./this->n_sides());
      bounds.second = 180. - (180./this->n_sides());
      break;

    case JACOBIAN:
    case SCALED_JACOBIAN:
      bounds.first  = 0.5;
      bounds.second = 1.;
      break;

    default:
      libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;

      bounds.second = -1;
    }

  return bounds;
}

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

      // Set any relevant node/edge maps for this element
      const std::string type_str (exio_helper->get_elem_type());
      const auto & conv = exio_helper->get_conversion(type_str);
      const bool is_c0polygon = (conv.libmesh_elem_type() == C0POLYGON);
      const bool is_c0polyhedron = (conv.libmesh_elem_type() == C0POLYHEDRON);
      std::size_t c0polygon_connect_offset = 0;
      std::size_t c0polyhedron_connect_offset = 0;

      // Loop over all the faces in this block
      int jmax = nelem_last_block+exio_helper->num_elem_this_blk;
      for (int j=nelem_last_block; j<jmax; j++)
        {
          const int elem_num = j - nelem_last_block;
          std::unique_ptr<Elem> uelem;

          if (is_c0polygon)
            {
              const int n_elem_nodes = exio_helper->elem_node_counts[elem_num];
              libmesh_error_msg_if(n_elem_nodes < 3,
                                   "Error: Exodus NSIDED block element "
                                   << elem_num
                                   << " has only " << n_elem_nodes << " nodes.");

              uelem = std::make_unique<C0Polygon>(cast_int<unsigned int>(n_elem_nodes));
            }
          else if (is_c0polyhedron)
            {
              const int n_faces = exio_helper->elem_face_counts[elem_num];
              libmesh_error_msg_if(n_faces < 4,
                                   "Error: Exodus NFACED block element "
                                   << elem_num
                                   << " has only " << n_faces << " faces.");

              std::vector<std::shared_ptr<Polygon>>
                sides(cast_int<std::size_t>(n_faces));

              for (auto s : index_range(sides))
                {
                  libmesh_assert_less(c0polyhedron_connect_offset,
                                      exio_helper->connect.size());
                  const int exodus_face_id =
                    exio_helper->connect[c0polyhedron_connect_offset++];

                  libmesh_error_msg_if
                    (exodus_face_id <= 0 ||
                     exodus_face_id >
                     cast_int<int>(exio_helper->c0polyhedron_face_connect.size()),

                     << " faces.");

                  const auto & face_nodes =
                    exio_helper->c0polyhedron_face_connect[exodus_face_id - 1];
                  auto side = std::make_shared<C0Polygon>
                    (cast_int<unsigned int>(face_nodes.size()));

                  for (auto n : index_range(face_nodes))
                    {
                      const auto libmesh_node_id =
                        exio_helper->get_libmesh_node_id(face_nodes[n]);
                      side->set_node(n, mesh.node_ptr(libmesh_node_id));
                    }

                  sides[s] = std::move(side);
                }

              std::unique_ptr<Node> mid_elem_node;
              uelem = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
              if (mid_elem_node)
                {
                  Node * added_node = mesh.add_node(std::move(mid_elem_node));
                  if (added_node->id() >= n_nodes)
                    n_nodes = added_node->id() + 1;
                }
            }
          else
            uelem = Elem::build(conv.libmesh_elem_type());

          // Make sure that Exodus's number of nodes per Elem matches
          // the number of Nodes for this type of Elem. We only check
          // this for the first Elem in each block, since these values
          // are the same for every Elem in the block.
          if (!is_c0polygon && !is_c0polyhedron && !elem_num)
            libmesh_error_msg_if(exio_helper->num_nodes_per_elem != static_cast<int>(uelem->n_nodes()),
                                 "Error: Exodus file says "
                                 << exio_helper->num_nodes_per_elem

          // If we don't have any Bezier extraction operators, this
          // is easy: we've already built all our nodes and just need
          // to link to them.
          if (is_c0polyhedron)
            {
              // Node pointers were already assigned while constructing
              // the polygonal sides above.
            }
          else if (exio_helper->bex_cv_conn.empty())
            {
              const int n_nodes_this_elem =
                is_c0polygon ?
                exio_helper->elem_node_counts[elem_num] :
                exio_helper->num_nodes_per_elem;

              for (int k=0; k<n_nodes_this_elem; k++)
                {
                  // Get index into this block's connectivity array
                  std::size_t gi =
                    is_c0polygon ?
                    c0polygon_connect_offset++ :
                    elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k);
                  libmesh_assert_less(gi, exio_helper->connect.size());

                  // Get the 1-based Exodus node id from the "connect" array
                  auto exodus_node_id = exio_helper->connect[gi];

            }
        }

      libmesh_assert(!is_c0polygon ||
                     c0polygon_connect_offset == exio_helper->connect.size());
      libmesh_assert(!is_c0polyhedron ||
                     c0polyhedron_connect_offset == exio_helper->connect.size());

      // running sum of # of elements per block,

  std::map<subdomain_id_type, std::vector<unsigned int>>
  build_subdomain_map(const MeshBase & mesh,
                      bool add_sides,
                      subdomain_id_type & subdomain_id_end,
                      int & next_block_id)

    // Allocate optional block IDs after both mesh subdomains and any blocks
    // synthesized for visualization sides.
    if (!subdomain_map.empty())
      next_block_id = cast_int<int>(subdomain_map.rbegin()->first) + 1;

    return subdomain_map;
  }

  num_elem_blk(header_info.num_elem_blk),
  num_edge(header_info.num_edge),
  num_edge_blk(header_info.num_edge_blk),
  num_face(header_info.num_face),
  num_face_blk(header_info.num_face_blk),
  num_node_sets(header_info.num_node_sets),
  num_side_sets(header_info.num_side_sets),
  num_elem_sets(header_info.num_elem_sets),

  convert_type(QUAD4, "QUAD4");
  convert_type(QUAD8, "QUAD8");
  convert_type(QUAD9, "QUAD9");
  convert_type(C0POLYGON, "NSIDED");
  {
    auto & conv = conversion_map[3][C0POLYHEDRON];
    conv.libmesh_type = C0POLYHEDRON;
    conv.exodus_type = "NFACED";
    conv.dim = 3;
    conv.n_nodes = 0;
  }
  convert_type(QUADSHELL4, "SHELL4", nullptr, nullptr, nullptr,
               /* inverse_side_map = */ &quadshell4_inverse_edge_map,

  element_equivalence_map["SHELL9"] = QUADSHELL9;

  // Runtime-topology polytope equivalences
  element_equivalence_map["NSIDED"] = C0POLYGON;
  element_equivalence_map["NFACED"] = C0POLYHEDRON;

  // TRI3 equivalences
  element_equivalence_map["TRI"]       = TRI3;

  if (read_only)
  {
    opened_for_reading = true;
    elem_node_counts.clear();
    elem_face_counts.clear();
    c0polyhedron_face_connect.clear();

    // ExodusII reads truncate to 32-char strings by default; we'd
    // like to support whatever's in the file, so as early as possible

  h.num_elem_sets = params.num_elem_sets;
  h.num_edge_blk = params.num_edge_blk;
  h.num_edge = params.num_edge;
  h.num_face_blk = params.num_face_blk;
  h.num_face = params.num_face;

  // And return it
  return h;

}


void ExodusII_IO_Helper::read_face_blocks()
{
  LOG_SCOPE("read_face_blocks()", "ExodusII_IO_Helper");

  if (!c0polyhedron_face_connect.empty())
    return;

  libmesh_error_msg_if(num_face_blk == 0,
                       "Error: Exodus NFACED element block found, "
                       "but the file has no face blocks.");

  std::vector<int> face_block_ids(num_face_blk);
  ex_err = exII::ex_get_ids(ex_id,
                            exII::EX_FACE_BLOCK,
                            face_block_ids.data());
  EX_CHECK_ERR(ex_err, "Error getting face block IDs.");

  c0polyhedron_face_connect.clear();
  c0polyhedron_face_connect.reserve(num_face);

  for (auto block : index_range(face_block_ids))
    {
      std::vector<char> face_type(libmesh_max_str_length+1);
      int num_face_this_blk = 0;
      int num_node_data_this_blk = 0;
      int num_edges_per_face = 0;
      int num_faces_per_face = 0;
      int num_attr_face = 0;

      ex_err = exII::ex_get_block(ex_id,
                                  exII::EX_FACE_BLOCK,
                                  face_block_ids[block],
                                  face_type.data(),
                                  &num_face_this_blk,
                                  &num_node_data_this_blk,
                                  &num_edges_per_face,
                                  &num_faces_per_face,
                                  &num_attr_face);
      EX_CHECK_ERR(ex_err, "Error getting face block info.");

      const auto & conv = get_conversion(std::string(face_type.data()));
      libmesh_error_msg_if(conv.libmesh_elem_type() != C0POLYGON,
                           "Error: NFACED polyhedron input currently expects "
                           "NSIDED face blocks, but face block "
                           << face_block_ids[block] << " has Exodus type "
                           << face_type.data() << ".");

      libmesh_error_msg_if
        (!(num_edges_per_face == 0) && !(num_edges_per_face == -1),
         "Error: Exodus NSIDED face block "
         << face_block_ids[block]
         << " has edge connectivity, which NFACED polyhedron input "
         << "does not currently support.");
      libmesh_error_msg_if
        (!(num_faces_per_face == 0) && !(num_faces_per_face == -1),
         "Error: Exodus NSIDED face block "
         << face_block_ids[block]
         << " has face-in-face connectivity, which NFACED polyhedron "
         << "input does not currently support.");

      std::vector<int> face_node_counts(num_face_this_blk);
      if (!face_node_counts.empty())
        {
          ex_err = exII::ex_get_entity_count_per_polyhedra
            (ex_id,
             exII::EX_FACE_BLOCK,
             face_block_ids[block],
             face_node_counts.data());
          EX_CHECK_ERR(ex_err, "Error reading polyhedron face node counts");
        }

      int counted_nodes = 0;
      for (const auto count : face_node_counts)
        counted_nodes += count;

      libmesh_error_msg_if(counted_nodes != num_node_data_this_blk,
                           "Error: Exodus NSIDED face block "
                           << face_block_ids[block]
                           << " says it has " << num_node_data_this_blk
                           << " total node entries, but its per-face "
                           << "node counts sum to " << counted_nodes << ".");

      std::vector<int> face_connect(num_node_data_this_blk);
      if (!face_connect.empty())
        {
          ex_err = exII::ex_get_conn(ex_id,
                                     exII::EX_FACE_BLOCK,
                                     face_block_ids[block],
                                     face_connect.data(),
                                     nullptr,
                                     nullptr);
          EX_CHECK_ERR(ex_err, "Error reading polyhedron face connectivity.");
        }

      std::size_t offset = 0;
      for (const auto count : face_node_counts)
        {
          libmesh_error_msg_if(count < 3,
                               "Error: Exodus NSIDED face block "
                               << face_block_ids[block]
                               << " has a face with only "
                               << count << " nodes.");

          c0polyhedron_face_connect.emplace_back
            (face_connect.begin() + offset,
             face_connect.begin() + offset + count);
          offset += count;
        }
    }

  libmesh_error_msg_if(c0polyhedron_face_connect.size() !=
                       cast_int<std::size_t>(num_face),
                       "Error: Exodus file says it has "
                       << num_face << " faces, but its face blocks contain "

  LOG_SCOPE("read_elem_in_block()", "ExodusII_IO_Helper");

  libmesh_assert_less (block, block_ids.size());
  elem_node_counts.clear();
  elem_face_counts.clear();

  // Unlike the other "extended" APIs, this one does not use a parameter struct.
  int num_edges_per_elem = 0;

  // Nemesis uses "Empty" as the element type for blocks with no
  // elements on the current processor.  There is no element conversion
  // for that sentinel type, nor is one needed for an empty block.
  const bool is_bezier = is_bezier_elem(elem_type.data());
  const Conversion * conversion = nullptr;
  if (num_elem_this_blk || is_bezier)
    conversion = &get_conversion(std::string(elem_type.data()));

  const bool is_c0polygon =
    conversion && conversion->libmesh_elem_type() == C0POLYGON;
  const bool is_c0polyhedron =
    conversion && conversion->libmesh_elem_type() == C0POLYHEDRON;

  // Warn or error when we don't currently support reading blocks with extended info.
  // Note: the docs say -1 will be returned for this but I found that it was
  // actually 0, so not sure which it will be in general.
  if (is_c0polyhedron && !(num_edges_per_elem == 0) && !(num_edges_per_elem == -1))
    libmesh_error_msg("Error: Exodus NFACED element blocks with edge "
                      "connectivity are not currently supported.");
  else if (!(num_edges_per_elem == 0) && !(num_edges_per_elem == -1))
    libmesh_warning("Exodus files with extended edge connectivity not currently supported.");
  if (!is_c0polyhedron && !(num_faces_per_elem == 0) && !(num_faces_per_elem == -1))
    libmesh_warning("Exodus files with extended face connectivity not currently supported.");

  // If we have a Bezier element here, then we've packed constraint

  // num_nodes_per_elem reflected both.
  if (is_bezier)
    {
      libmesh_assert(conversion);
      num_nodes_per_elem = conversion->n_nodes;
    }
  else if (is_c0polygon)
    {
      elem_node_counts.resize(num_elem_this_blk);

      if (!elem_node_counts.empty())
        {
          ex_err = exII::ex_get_entity_count_per_polyhedra
            (ex_id,
             exII::EX_ELEM_BLOCK,
             block_ids[block],
             elem_node_counts.data());
          EX_CHECK_ERR(ex_err, "Error reading polygon node counts");
        }

      int counted_nodes = 0;
      for (const auto count : elem_node_counts)
        counted_nodes += count;

      libmesh_error_msg_if(counted_nodes != num_node_data_per_elem,
                           "Error: Exodus NSIDED block "
                           << block_ids[block]
                           << " says it has " << num_node_data_per_elem
                           << " total node entries, but its per-element "
                           << "node counts sum to " << counted_nodes << ".");

      num_nodes_per_elem = 0;
    }
  else if (is_c0polyhedron)
    {
      if (c0polyhedron_face_connect.empty())
        this->read_face_blocks();

      elem_face_counts.resize(num_elem_this_blk);

      if (!elem_face_counts.empty())
        {
          ex_err = exII::ex_get_entity_count_per_polyhedra
            (ex_id,
             exII::EX_ELEM_BLOCK,
             block_ids[block],
             elem_face_counts.data());
          EX_CHECK_ERR(ex_err, "Error reading polyhedron face counts");
        }

      int counted_faces = 0;
      for (const auto count : elem_face_counts)
        counted_faces += count;

      libmesh_error_msg_if(counted_faces != num_faces_per_elem,
                           "Error: Exodus NFACED block "
                           << block_ids[block]
                           << " says it has " << num_faces_per_elem
                           << " total face entries, but its per-element "
                           << "face counts sum to " << counted_faces << ".");

      num_nodes_per_elem = 0;
    }
  else
    num_nodes_per_elem = num_node_data_per_elem;

  if (verbose)
    {
      libMesh::out << "Read a block of " << num_elem_this_blk
                   << " " << elem_type.data() << "(s)";
      if (is_c0polygon)
        libMesh::out << " having " << num_node_data_per_elem
                     << " total node entries.";
      else if (is_c0polyhedron)
        libMesh::out << " having " << num_faces_per_elem
                     << " total face entries.";
      else
        libMesh::out << " having " << num_nodes_per_elem
                     << " nodes per element.";
      libMesh::out << std::endl;
    }

  // Read in the connectivity of the elements of this block,
  // watching out for the case where we actually have no
  // elements in this block (possible with parallel files)
  connect.resize(is_c0polygon ?
                 num_node_data_per_elem :
                 is_c0polyhedron ?
                 num_faces_per_elem :
                 num_node_data_per_elem*num_elem_this_blk);

  if (!connect.empty())
    {
      ex_err = exII::ex_get_conn(ex_id,
                                 exII::EX_ELEM_BLOCK,
                                 block_ids[block],
                                 is_c0polyhedron ? nullptr : connect.data(),
                                 nullptr, // elem_edge_conn (unused)
                                 is_c0polyhedron ? connect.data() : nullptr);

      EX_CHECK_ERR(ex_err, "Error reading block connectivity.");
      message("Connectivity retrieved successfully for block: ", block);

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

  // We need to know about all processors' subdomains
  subdomain_id_type subdomain_id_end = 0;
  int c0polyhedron_face_block_id = -1;
  auto subdomain_map = build_subdomain_map(mesh,
                                           _add_sides,
                                           subdomain_id_end,
                                           c0polyhedron_face_block_id);

  num_elem = n_active_elem;
  num_nodes = 0;
  num_face = 0;
  num_face_blk = 0;

  dof_id_type local_num_c0polyhedron_faces = 0;
  bool has_c0polyhedron = false;
  for (const auto & elem : mesh.active_local_element_ptr_range())
    if (elem->type() == C0POLYHEDRON)
      {
        has_c0polyhedron = true;
        local_num_c0polyhedron_faces += elem->n_sides();
      }

  mesh.comm().sum(local_num_c0polyhedron_faces);
  mesh.comm().max(has_c0polyhedron);
  if (has_c0polyhedron)
    {
      num_face = cast_int<int>(local_num_c0polyhedron_faces);
      num_face_blk = 1;
    }

  // If we're adding face elements they'll need copies of their nodes.

  params.num_elem_sets = num_elem_sets;
  params.num_edge_blk = num_edge_blk;
  params.num_edge = num_edge;
  params.num_face_blk = num_face_blk;
  params.num_face = num_face;

  ex_err = exII::ex_put_init_ext(ex_id, &params);
  EX_CHECK_ERR(ex_err, "Error initializing new Exodus file.");

  // Map from block ID to a vector of element IDs in that block.  Element
  // IDs are now of type dof_id_type, subdomain IDs are of type subdomain_id_type.
  subdomain_id_type subdomain_id_end = 0;
  int c0polyhedron_face_block_id = -1;
  auto subdomain_map = build_subdomain_map(mesh,
                                           _add_sides,
                                           subdomain_id_end,
                                           c0polyhedron_face_block_id);

  if ((_run_only_on_proc0) && (this->processor_id() != 0))
    return;

  NamesData names_table(num_elem_blk, _max_name_length);

  num_elem = 0;
  bool has_c0polygon_blocks = false;
  bool has_c0polyhedron_blocks = false;
  int c0polyhedron_total_faces = 0;
  int c0polyhedron_total_face_nodes = 0;

  // counter indexes into the block_ids vector
  unsigned int counter = 0;

          num_elem_this_blk_vec.push_back
            (cast_int<int>(element_id_vec.size()));

          std::string block_name = mesh.subdomain_name(subdomain_id);
          if (block_name.empty() && elem_t == C0POLYGON)
            block_name = "NSIDED_" + std::to_string(counter + 1);
          if (block_name.empty() && elem_t == C0POLYHEDRON)
            block_name = "NFACED_" + std::to_string(counter + 1);
          names_table.push_back_entry(block_name);
        }

      num_elem += num_elem_this_blk_vec.back();

      // Note that Exodus assumes all elements in a block are of the same type!
      // We are using that same assumption here!
      const auto & conv = get_conversion(elem_t);
      int num_edges_per_elem = 0;
      int num_faces_per_elem = 0;
      if (elem_t == C0POLYGON)
        {
          if (subdomain_id >= subdomain_id_end)
            libmesh_not_implemented_msg("Support for C0POLYGON side blocks not yet implemented");

          has_c0polygon_blocks = true;
          num_nodes_per_elem = 0;

          for (auto elem_id : element_id_vec)
            {
              const Elem & elem = mesh.elem_ref(elem_id);

              libmesh_error_msg_if(elem.type() != C0POLYGON,
                                   "Error: Exodus requires all elements with a given subdomain ID "
                                   "to be the same type.\n"
                                   << "Can't write both "
                                   << Utility::enum_to_string(elem.type())
                                   << " and C0POLYGON in the same block!");

              num_nodes_per_elem += cast_int<int>(elem.n_nodes());
            }
        }
      else if (elem_t == C0POLYHEDRON)
        {
          if (subdomain_id >= subdomain_id_end)
            libmesh_not_implemented_msg("Support for C0POLYHEDRON side blocks not yet implemented");

          has_c0polyhedron_blocks = true;
          num_nodes_per_elem = 0;

          for (auto elem_id : element_id_vec)
            {
              const Elem & elem = mesh.elem_ref(elem_id);

              libmesh_error_msg_if(elem.type() != C0POLYHEDRON,
                                   "Error: Exodus requires all elements with a given subdomain ID "
                                   "to be the same type.\n"
                                   << "Can't write both "
                                   << Utility::enum_to_string(elem.type())
                                   << " and C0POLYHEDRON in the same block!");

              const int elem_n_sides = cast_int<int>(elem.n_sides());
              num_faces_per_elem += elem_n_sides;
              c0polyhedron_total_faces += elem_n_sides;

              for (auto s : elem.side_index_range())
                c0polyhedron_total_face_nodes += cast_int<int>(elem.nodes_on_side(s).size());
            }
        }
      else
        {
          num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t];
          if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint)
            libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
        }

      elem_blk_id.push_back(subdomain_id);
      elem_type_table.push_back_entry(conv.exodus_elem_type().c_str());
      num_nodes_per_elem_vec.push_back(num_nodes_per_elem);
      num_attr_vec.push_back(0); // we don't currently use elem block attributes.
      num_edges_per_elem_vec.push_back(num_edges_per_elem); // We don't currently store any edge blocks
      num_faces_per_elem_vec.push_back(num_faces_per_elem);
      ++counter;
    }

  if (has_c0polyhedron_blocks)
    {
      libmesh_assert_equal_to(num_face_blk, 1);
      libmesh_assert_equal_to(num_face, c0polyhedron_total_faces);
    }

  // Here we reserve() space so that we can push_back() onto the

  // We also build a data structure of edge block names which can
  // later be passed to exII::ex_put_names().
  NamesData edge_block_names_table(num_edge_blk, _max_name_length);
  NamesData face_block_names_table(num_face_blk, _max_name_length);
  if (has_c0polyhedron_blocks)
    face_block_names_table.push_back_entry("NSIDED_FACES");

  // Note: We are going to use the edge **boundary** ids as **block** ids.
  for (const auto & pr : edge_id_to_conn)

      edge_block_names_table.push_back_entry(bi.get_edgeset_name(id));
    }

  if (has_c0polygon_blocks || has_c0polyhedron_blocks)
    {
      // ex_put_concat_all_blocks() does not define the per-polytope
      // entity count arrays required by NSIDED/NFACED blocks in all supported
      // Exodus versions. Define blocks individually through ex_put_block().
      if (has_c0polyhedron_blocks)
        {
          ex_err = exII::ex_put_block(ex_id,
                                      exII::EX_FACE_BLOCK,
                                      c0polyhedron_face_block_id,
                                      "NSIDED",

                                      0,
                                      0,
                                      0);
          EX_CHECK_ERR(ex_err, "Error writing polyhedron face block.");
        }

      for (auto i : index_range(elem_blk_id))
        {
          ex_err = exII::ex_put_block(ex_id,
                                      exII::EX_ELEM_BLOCK,
                                      elem_blk_id[i],
                                      elem_type_table.get_char_star(cast_int<int>(i)),
                                      num_elem_this_blk_vec[i],
                                      num_nodes_per_elem_vec[i],
                                      num_edges_per_elem_vec[i],
                                      num_faces_per_elem_vec[i],
                                      num_attr_vec[i]);
          EX_CHECK_ERR(ex_err, "Error writing element block.");
        }

      for (auto i : index_range(edge_blk_id))
        {
          ex_err = exII::ex_put_block(ex_id,
                                      exII::EX_EDGE_BLOCK,
                                      edge_blk_id[i],
                                      edge_type_table.get_char_star(cast_int<int>(i)),
                                      num_edge_this_blk_vec[i],
                                      num_nodes_per_edge_vec[i],
                                      0,
                                      0,
                                      num_attr_edge_vec[i]);
          EX_CHECK_ERR(ex_err, "Error writing edge block.");
        }
    }
  else
    {

      // problem is that we need to first specify a nonzero number of
      // node/elem maps during the call to ex_put_init_ext() in order for
      // this to work correctly.
      exII::ex_block_params params = {};

      // Set pointers for information about elem blocks.
      params.elem_blk_id = elem_blk_id.data();
      params.elem_type = elem_type_table.get_char_star_star();
      params.num_elem_this_blk = num_elem_this_blk_vec.data();
      params.num_nodes_per_elem = num_nodes_per_elem_vec.data();
      params.num_edges_per_elem = num_edges_per_elem_vec.data();
      params.num_faces_per_elem = num_faces_per_elem_vec.data();
      params.num_attr_elem = num_attr_vec.data();
      params.define_maps = 0;

      // Set pointers to edge block information only if we actually have some.
      if (num_edge_blk)
        {
          params.edge_blk_id = edge_blk_id.data();
          params.edge_type = edge_type_table.get_char_star_star();
          params.num_edge_this_blk = num_edge_this_blk_vec.data();
          params.num_nodes_per_edge = num_nodes_per_edge_vec.data();
          params.num_attr_edge = num_attr_edge_vec.data();
        }

      ex_err = exII::ex_put_concat_all_blocks(ex_id, &params);
      EX_CHECK_ERR(ex_err, "Error writing element blocks.");
    }

  // This counter is used to fill up the libmesh_elem_num_to_exodus map in the loop below.

    next_fake_id = mesh.next_unique_id();
#endif

  std::vector<int> face_connect;
  std::vector<int> c0polyhedron_face_node_counts;
  face_connect.reserve(c0polyhedron_total_face_nodes);
  c0polyhedron_face_node_counts.reserve(c0polyhedron_total_faces);
  int next_c0polyhedron_face_id = 1;

  const auto get_exodus_node_id = [&](const Elem &elem,
                                  dof_id_type elem_id,
                                  unsigned int elem_node_index)
                                   -> int
  {
    if (!use_discontinuous)
    return libmesh_map_find(libmesh_node_num_to_exodus,
                            cast_int<int>(elem.node_id(elem_node_index)));

    return cast_int<int>(libmesh_map_find(discontinuous_node_indices,
                                          std::make_pair(elem_id, elem_node_index)));
  };

  for (auto & [subdomain_id, element_id_vec] : subdomain_map)
    {

        mesh.elem_ref(element_id_vec[0]).type();

      const auto & conv = get_conversion(elem_t);
      const bool is_c0polygon_block = (elem_t == C0POLYGON);
      const bool is_c0polyhedron_block = (elem_t == C0POLYHEDRON);
      const bool is_variable_connectivity_block =
        is_c0polygon_block || is_c0polyhedron_block;
      std::vector<int> c0polygon_node_counts;
      std::vector<int> c0polyhedron_face_counts;

      if (is_variable_connectivity_block && subdomain_id >= subdomain_id_end)
        libmesh_not_implemented_msg("Support for " <<
                                    Utility::enum_to_string(elem_t) <<
                                    " side blocks not yet implemented");

      if (!is_variable_connectivity_block)
        {
          num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t];
          if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint)
            libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented");
        }

      // If this is a *real* block, we just loop over vectors of
      // element ids to add.
      if (subdomain_id < subdomain_id_end)
      {
        if (is_variable_connectivity_block)
            connect.clear();
        if (is_c0polygon_block)
            c0polygon_node_counts.reserve(element_id_vec.size());
        else if (is_c0polyhedron_block)
            c0polyhedron_face_counts.reserve(element_id_vec.size());
        else
          connect.resize(element_id_vec.size()*num_nodes_per_elem);

        const auto add_c0polygon_connectivity =
            [&](const Elem &elem, dof_id_type elem_id)
        {
          c0polygon_node_counts.push_back(cast_int<int>(elem.n_nodes()));
          for (auto elem_node_index : elem.node_index_range())
            connect.push_back(get_exodus_node_id(elem, elem_id, elem_node_index));

        };

        const auto add_c0polyhedron_connectivity =
            [&](const Elem &elem, dof_id_type elem_id)
        {
          c0polyhedron_face_counts.push_back(cast_int<int>(elem.n_sides()));

          for (auto s: elem.side_index_range())
          {
            connect.push_back(next_c0polyhedron_face_id++);

            const std::vector<unsigned int> side_nodes = elem.nodes_on_side(s);
            c0polyhedron_face_node_counts.push_back(cast_int<int>(side_nodes.size()));

            for (const auto elem_node_index : side_nodes)
            face_connect.push_back(
              get_exodus_node_id(elem, elem_id, elem_node_index));
          }
        };

        const auto add_fixed_connectivity =
          [&](const Elem &elem, dof_id_type elem_id, std::size_t elem_index)
        {
          for (unsigned int j = 0;
              j < static_cast<unsigned int>(num_nodes_per_elem);
              ++j)
          {
            const auto connect_index =
              cast_int<unsigned int>((elem_index * num_nodes_per_elem) + j);
            const auto elem_node_index = conv.get_inverse_node_map(j);
            connect[connect_index] =
              get_exodus_node_id(elem, elem_id, elem_node_index);
          }
        };

        for (auto i : index_range(element_id_vec))
        {

                               << Utility::enum_to_string(conv.libmesh_elem_type())
                               << " in the same block!");

          if (is_c0polygon_block)
            add_c0polygon_connectivity(elem, elem_id);
          else if (is_c0polyhedron_block)
            add_c0polyhedron_connectivity(elem, elem_id);
          else
            add_fixed_connectivity(elem, elem_id, i);

          // push_back() either elem_id+1 or the current Elem's
          // unique_id+1 into the elem_num_map, depending on the value

        (ex_id,
         exII::EX_ELEM_BLOCK,
         subdomain_id,
         is_c0polyhedron_block ? nullptr : connect.data(), // node_conn
         nullptr,                                          // elem_edge_conn (unused)
         is_c0polyhedron_block ? connect.data() : nullptr);
      EX_CHECK_ERR(ex_err, "Error writing element connectivities");

      if (is_c0polygon_block)
        {
          ex_err = exII::ex_put_entity_count_per_polyhedra
            (ex_id,
             exII::EX_ELEM_BLOCK,
             subdomain_id,
             c0polygon_node_counts.data());
          EX_CHECK_ERR(ex_err, "Error writing polygon node counts");
        }
      if (is_c0polyhedron_block)
        {
          ex_err = exII::ex_put_entity_count_per_polyhedra
            (ex_id,
             exII::EX_ELEM_BLOCK,
             subdomain_id,
             c0polyhedron_face_counts.data());
          EX_CHECK_ERR(ex_err, "Error writing polyhedron face counts");
        }
    } // end for (auto & [subdomain_id, element_id_vec] : subdomain_map)

  if (has_c0polyhedron_blocks)
    {
      libmesh_assert_equal_to(c0polyhedron_face_node_counts.size(),
                              cast_int<std::size_t>(num_face));
      libmesh_assert_equal_to(next_c0polyhedron_face_id, num_face + 1);

      ex_err = exII::ex_put_conn
        (ex_id,
         exII::EX_FACE_BLOCK,
         c0polyhedron_face_block_id,
         face_connect.data(), // node_conn
         nullptr,             // elem_edge_conn (unused)
         nullptr);            // elem_face_conn (unused)
      EX_CHECK_ERR(ex_err, "Error writing polyhedron face connectivities");

      ex_err = exII::ex_put_entity_count_per_polyhedra
        (ex_id,
         exII::EX_FACE_BLOCK,
         c0polyhedron_face_block_id,
         c0polyhedron_face_node_counts.data());
      EX_CHECK_ERR(ex_err, "Error writing polyhedron face node counts");
    }

  // write out the element number map that we created

      EX_CHECK_ERR(ex_err, "Error writing element block names");
    }

  if (num_face_blk > 0)
    {
      ex_err = exII::ex_put_names
        (ex_id,
         exII::EX_FACE_BLOCK,
         face_block_names_table.get_char_star_star());
      EX_CHECK_ERR(ex_err, "Error writing face block names");
    }

  // Write out edge blocks if we have any

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

            std::vector<Node *> node_list;

            // Start with a layer of triangles on the boundary
            const auto dx_tri = Real(1) / nx;
            const auto dy_tri = Real(1) / (ny + 1);
            std::unique_ptr<Elem> new_elem;
            for (const auto i : make_range(nx + 1))
            {

            // Build layers of hexagons
            const auto hex_side =
                (Real(1) - (ny == 1 ?
                             dy_tri :
                             (Real(1) + (ny - 1) / 2.) * dy_tri)) / ny;
            for (const auto j : make_range(ny))
            {
              for (const auto i : make_range(nx + (j % 2)))

              Node *node0, *node1, *node2;
              if (i == 0 && ny_odd)
              {
                node0 = mesh.add_point(Point(0., 1., 0.));
                node1 = node_list[running_index++];
                node2 = node_list[running_index];
              }

              {
                node0 = node_list[running_index++];
                node1 = node_list[running_index];
                node2 = mesh.add_point(Point(1., 1., 0.));
              }

              new_elem = std::make_unique<C0Polygon>(3);

          case C0POLYHEDRON:
            {
              const dof_id_type grid_nodes =
                cast_int<dof_id_type>((nx+1)*(ny+1)*(nz+1));

              // Reserve one interior node per polyhedron for the robust
              // fallback tetrahedralization used when the preferred
              // tetrahedralization cannot be constructed.
              const dof_id_type mid_polyhedron_nodes =
                (type == C0POLYHEDRON) ?
                cast_int<dof_id_type>(nx*ny*nz) : 0;

              mesh.reserve_nodes(grid_nodes + mid_polyhedron_nodes);
              break;
            }

                      elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
                      elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  ));

                      if (k == 0)
                        boundary_info.add_side(elem, 0, 0);

                      if (k == (nz-1))
                        boundary_info.add_side(elem, 5, 5);

                      if (j == 0)
                        boundary_info.add_side(elem, 1, 1);

                      if (j == (ny-1))
                        boundary_info.add_side(elem, 3, 3);

                      if (i == 0)
                        boundary_info.add_side(elem, 4, 4);

                      if (i == (nx-1))
                        boundary_info.add_side(elem, 2, 2);
                    }
              break;
            }


          case C0POLYHEDRON:
            {
              const std::array<std::array<unsigned int, 4>, 6> side_nodes =
                {{{0, 1, 2, 3},   // min z
                  {0, 1, 5, 4},   // min y
                  {2, 6, 5, 1},   // max x

                  {0, 4, 7, 3},   // min x
                  {5, 6, 7, 4}}}; // max z

              for (unsigned int k=0; k<nz; k++)
                for (unsigned int j=0; j<ny; j++)
                  for (unsigned int i=0; i<nx; i++)
                    {
                      std::array<Node *, 8> elem_nodes =
                        {{mesh.node_ptr(idx(type,nx,ny,i,j,k)      ),
                          mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    ),
                          mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  ),
                          mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    ),
                          mesh.node_ptr(idx(type,nx,ny,i,j,k+1)    ),
                          mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  ),
                          mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)),
                          mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  )}};

                      std::vector<std::shared_ptr<Polygon>> sides(side_nodes.size());
                      for (auto s : index_range(side_nodes))
                        {
                          sides[s] = std::make_shared<C0Polygon>(side_nodes[s].size());
                          for (auto n : index_range(side_nodes[s]))
                            sides[s]->set_node(n, elem_nodes[side_nodes[s][n]]);
                        }

                      std::unique_ptr<Node> mid_elem_node;
                      std::unique_ptr<Elem> new_elem =
                        std::make_unique<C0Polyhedron>(sides, mid_elem_node);
                      if (mid_elem_node)
                        mesh.add_node(std::move(mid_elem_node));

                      new_elem->set_id() = elem_id++;
                      Elem * elem = mesh.add_elem(std::move(new_elem));

                      if (k == 0)
                        boundary_info.add_side(elem, 0, 0);

                      if (i == (nx-1))
                        boundary_info.add_side(elem, 2, 2);
                    }
              break;
            }

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

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

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

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

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

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

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

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

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

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

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

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

          elems_to_query[pid].push_back(coarse_id);
        }

      // Return fine element data based on vertex_average()

        elem_refinement_datum;

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

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

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

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

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

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

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

}


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

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

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

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

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

          hi_node->processor_id() = chosen_pid;
        }

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

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

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

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

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

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

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

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

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

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

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

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

              bool child_contains_our_nodes = true;

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

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

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

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

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

    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

  const ElemType type =
    cast_int<ElemType>(typeint);

  unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
  unsigned int n_sides = Elem::type_to_n_sides_map[type];
  unsigned int n_edges = Elem::type_to_n_edges_map[type];
  unsigned int variable_topology_size = 0;

  if (n_nodes == invalid_uint)
    {
      libmesh_error_msg_if(type != C0POLYGON && type != C0POLYHEDRON,
                           "Unsupported variable-topology element "
                           << Utility::enum_to_string(type));

      auto topology_in = in + header_size;
      n_nodes = cast_int<unsigned int>(*topology_in++);
      n_sides = cast_int<unsigned int>(*topology_in++);
      n_edges = cast_int<unsigned int>(*topology_in++);
      variable_topology_size = 3;

      if (type == C0POLYHEDRON)
        {
          topology_in += n_nodes;
          for (unsigned int s = 0; s != n_sides; ++s)
            {
              const unsigned int n_side_nodes =
                cast_int<unsigned int>(*topology_in++);
              topology_in += n_side_nodes;
              variable_topology_size += n_side_nodes + 1;
            }
        }
    }

  const unsigned int pre_indexing_size =
    header_size + variable_topology_size + n_nodes + n_sides*2;

  const unsigned int indexing_size =
    DofObject::unpackable_indexing_size(in+pre_indexing_size);

Packing<const Elem *>::packable_size (const Elem * const & elem,
                                      const MeshBase * mesh)
{
  unsigned int variable_topology_size = 0;
  if (elem->type() == C0POLYGON || elem->type() == C0POLYHEDRON)
    {
      // Store the dynamic node, side, and edge counts.
      variable_topology_size = 3;

      // A polygon's node ordering fully specifies its topology.  A
      // polyhedron additionally needs each side's local node indices.
      if (elem->type() == C0POLYHEDRON)
        for (auto s : elem->side_index_range())
          variable_topology_size +=
            1 + cast_int<unsigned int>(elem->nodes_on_side(s).size());
    }

  // We always communicate if we are on a boundary or not
  unsigned int total_packed_bcs = 1;
  const unsigned int n_sides = elem->n_sides();

  largest_id_type on_boundary = 0;
  for (auto s : elem->side_index_range())

    if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side())
    {
      total_packed_bcs += n_sides;
      for (unsigned int s = 0; s != n_sides; ++s)
        total_packed_bcs +=
          mesh->get_boundary_info().n_raw_boundary_ids(elem,s);
    }

  if (elem->level() == 0)
    {
      const unsigned int n_edges = elem->n_edges();
      total_packed_bcs += n_edges;
      for (unsigned int e = 0; e != n_edges; ++e)
        total_packed_bcs +=
          mesh->get_boundary_info().n_edge_boundary_ids(elem,e);

#ifndef NDEBUG
    1 + // add an int for the magic header when testing
#endif
    header_size + variable_topology_size + elem->n_nodes() + n_sides*2 +
    elem->packed_indexing_size() + total_packed_bcs;
}

    *data_out++ =(DofObject::invalid_id);

  const bool has_variable_topology =
    elem->type() == C0POLYGON || elem->type() == C0POLYHEDRON;
  if (has_variable_topology)
    {
      *data_out++ = elem->n_nodes();
      *data_out++ = elem->n_sides();
      *data_out++ = elem->n_edges();
    }

  for (const Node & node : elem->node_ref_range())
    *data_out++ = node.id();

  if (elem->type() == C0POLYHEDRON)
    for (auto s : elem->side_index_range())
      {
        const std::vector<unsigned int> side_nodes =
          elem->nodes_on_side(s);
        *data_out++ = side_nodes.size();
        for (const auto node : side_nodes)
          *data_out++ = node;
      }

  // Add the id of and the side for any return link from each neighbor

  const ElemType type =
    cast_int<ElemType>(typeint);

  unsigned int n_nodes = Elem::type_to_n_nodes_map[type];
  unsigned int n_sides = Elem::type_to_n_sides_map[type];
  unsigned int n_edges = Elem::type_to_n_edges_map[type];

  // int 5: processor id
  const processor_id_type processor_id =

  // plus the real data header
  libmesh_assert_equal_to (in - original_in, header_size + 1);

  if (n_nodes == invalid_uint)
    {
      libmesh_error_msg_if(type != C0POLYGON && type != C0POLYHEDRON,
                           "Unsupported variable-topology element "
                           << Utility::enum_to_string(type));

      n_nodes = cast_int<unsigned int>(*in++);
      n_sides = cast_int<unsigned int>(*in++);
      n_edges = cast_int<unsigned int>(*in++);

      if (type == C0POLYGON)
        libmesh_error_msg_if
          (n_nodes < 3 || n_sides != n_nodes || n_edges != n_nodes,
           "Invalid packed C0POLYGON topology with "
           << n_nodes << " nodes, " << n_sides << " sides, and "
           << n_edges << " edges");
      else
        libmesh_error_msg_if
          (n_nodes < 4 || n_sides < 4,
           "Invalid packed C0POLYHEDRON topology with "
           << n_nodes << " nodes and " << n_sides << " sides");
    }

  const auto node_ids_in = in;
  in += n_nodes;

  std::vector<std::vector<unsigned int>> polyhedron_side_nodes;
  if (type == C0POLYHEDRON)
    {
      polyhedron_side_nodes.resize(n_sides);
      std::vector<bool> node_seen(n_nodes, false);
      unsigned int next_new_node = 0;
      for (auto & side_nodes : polyhedron_side_nodes)
        {
          const unsigned int n_side_nodes =
            cast_int<unsigned int>(*in++);
          libmesh_error_msg_if(n_side_nodes < 3,
                               "Cannot unpack a C0POLYHEDRON side with only "
                               << n_side_nodes << " nodes");
          side_nodes.resize(n_side_nodes);
          for (auto & node : side_nodes)
            {
              node = cast_int<unsigned int>(*in++);
              libmesh_error_msg_if(node >= n_nodes,
                                   "Invalid local node " << node
                                   << " in packed C0POLYHEDRON side");

              if (!node_seen[node])
                {
                  libmesh_error_msg_if
                    (node != next_new_node,
                     "Packed C0POLYHEDRON side connectivity first encounters "
                     "local node " << node << " where local node "
                     << next_new_node << " was expected");
                  node_seen[node] = true;
                  ++next_new_node;
                }
            }
        }

      libmesh_error_msg_if
        (next_new_node != n_nodes && next_new_node + 1 != n_nodes,
         "Packed C0POLYHEDRON sides reference " << next_new_node
         << " vertices, but the element has " << n_nodes

      libmesh_assert_equal_to (elem->subdomain_id(), subdomain_id);
      libmesh_assert_equal_to (elem->type(), type);

      if (type == C0POLYGON || type == C0POLYHEDRON)
        {
          libmesh_error_msg_if
            (elem->type() != type ||
             elem->n_nodes() != n_nodes ||
             elem->n_sides() != n_sides ||

             "Existing " << Utility::enum_to_string(type)
             << " topology does not match its packed topology");

          for (unsigned int n = 0; n != n_nodes; ++n)
            libmesh_error_msg_if
              (elem->node_id(n) !=
               cast_int<dof_id_type>(*(node_ids_in + n)),
               "Existing " << Utility::enum_to_string(type)
               << " local node " << n
               << " does not match its packed node");

          if (type == C0POLYHEDRON)
            for (auto s : elem->side_index_range())
              libmesh_error_msg_if
                (elem->nodes_on_side(s) != polyhedron_side_nodes[s],
                 "Existing C0POLYHEDRON side " << s
                 << " does not match its packed topology");
        }

      libmesh_assert_equal_to (elem->n_nodes(), n_nodes);
      libmesh_assert_equal_to (elem->n_sides(), n_sides);
      libmesh_assert_equal_to (elem->n_edges(), n_edges);

#ifndef NDEBUG
      // All our nodes should be correct

      libmesh_assert_equal_to (level, 0);
#endif

      if (type == C0POLYGON)
        elem = std::make_unique<C0Polygon>(n_nodes, parent).release();
      else if (type == C0POLYHEDRON)
        {
          std::vector<std::shared_ptr<Polygon>> sides(n_sides);
          for (auto s : index_range(sides))
            {
              const auto & side_nodes = polyhedron_side_nodes[s];
              auto side = std::make_shared<C0Polygon>
                (cast_int<unsigned int>(side_nodes.size()));
              for (auto n : index_range(side_nodes))
                {
                  const dof_id_type node_id =
                    cast_int<dof_id_type>
                    (*(node_ids_in + side_nodes[n]));
                  side->set_node(n, mesh->node_ptr(node_id));
                }
              sides[s] = std::move(side);
            }

          std::unique_ptr<Node> generated_mid_node;
          auto polyhedron = std::make_unique<C0Polyhedron>
            (sides, generated_mid_node, parent);

          libmesh_error_msg_if
            (polyhedron->n_nodes() != n_nodes,
             "Packed C0POLYHEDRON has " << n_nodes
             << " nodes, but reconstructing its topology produced "
             << polyhedron->n_nodes() << " nodes");

          if (generated_mid_node)
            {
              const dof_id_type mid_node_id =
                cast_int<dof_id_type>(*(node_ids_in + n_nodes - 1));
              polyhedron->set_node(n_nodes - 1,
                                   mesh->node_ptr(mid_node_id));
            }

          elem = polyhedron.release();
        }
      else
        elem = Elem::build(type,parent).release();
      libmesh_assert (elem);

#ifdef LIBMESH_ENABLE_AMR

      // Assign the connectivity
      libmesh_assert_equal_to (elem->n_nodes(), n_nodes);

      if (type == C0POLYHEDRON)
        {
          libmesh_error_msg_if(elem->n_sides() != n_sides,
                               "Packed C0POLYHEDRON has " << n_sides
                               << " sides, but reconstructing its topology produced "
                               << elem->n_sides() << " sides");
          libmesh_error_msg_if(elem->n_edges() != n_edges,
                               "Packed C0POLYHEDRON has " << n_edges
                               << " edges, but reconstructing its topology produced "
                               << elem->n_edges() << " edges");

          for (unsigned int n = 0; n != n_nodes; ++n)
            libmesh_error_msg_if
              (elem->node_id(n) !=
               cast_int<dof_id_type>(*(node_ids_in + n)),
               "Packed C0POLYHEDRON local node ordering was not preserved");
        }
      else
        for (unsigned int n=0; n != n_nodes; n++)
          elem->set_node
            (n, mesh->node_ptr
             (cast_int<dof_id_type>(*(node_ids_in + n))));

      // Set interior_parent if found
      {

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

  // Return the new element
  return elem;
}

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