* values of the finite element solution at the vertices, and scaling them would change what a nodal
 * boundary condition or a nodal output of a HIERARCHIC variable means.
 */
inline Real fe_hierarchic_bubble_scaling(const unsigned int i)
{
  libmesh_assert_greater(i, 1);

  // An even bubble differentiates to xi^(i-1)/(i-1)!, whose square integrates to 2/(2i-1) over the
  // interval. The linear term an odd bubble carries turns the i^2 of that calculation into (i-1)^2.
  const Real denominator = (i % 2) ? Real(i) - 1. : Real(i);

  return std::sqrt((2. * Real(i) - 1.) / 2.) / denominator;
}

  /**
   * \returns 5.  All pyramids have 5 faces.
   */
  virtual unsigned int n_faces() const override { return 5; }

  /**
   * \returns 10.

  PetscMatrixBase & operator= (PetscMatrixBase &&) = delete;
  virtual ~PetscMatrixBase ();

  virtual SolverPackage solver_package() override
  {
    return PETSC_SOLVERS;
  }

  /**

template <typename T>
void
PetscMFFDMatrix<T>::assign(Mat m, bool set_context)
{
  this->_mat = m;
  this->_is_initialized = true;
  this->_destroy_mat_on_exit = false;
  if (set_context)
    this->set_context();
}

template <typename T>
void

   * the preconditioning matrix -- which may be the same object (the common case) or genuinely
   * distinct (e.g. a matrix-free operator paired with an assembled preconditioning matrix).
   */
  virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> & /* jac_in */,
                                               SparseMatrix<T> & /* pre_in */,
                                               NumericVector<T> & /* x_in */,
                                               NumericVector<T> & /* r_in */,
                                               const double /* tol */,
                                               const unsigned int /* m_its */)
  {
    libmesh_not_implemented();
  }

  /**

// T::operator*) might do worse, so we'll specialize here.
template <typename T>
struct do_pow<6,T> {
  static inline T apply (const T & x)
  {
    const T x2 = x*x,
      x4 = x2*x2;

    return x4*x2;
  }
};

class AddConstraint
{
public:
  virtual ~AddConstraint() = default;

  virtual void operator()(dof_id_type dof_number,
                          const DofConstraintRow & constraint_row,
                          const Number constraint_rhs) const = 0;
};

class AddPrimalConstraint : public AddConstraint
{
private:
  DofMap                  & dof_map;

public:
  AddPrimalConstraint(DofMap & dof_map_in) : dof_map(dof_map_in) {}

  virtual void operator()(dof_id_type dof_number,
                          const DofConstraintRow & constraint_row,

public:
  AddAdjointConstraint(DofMap & dof_map_in, unsigned int qoi_index_in)
    : dof_map(dof_map_in), qoi_index(qoi_index_in) {}

  virtual void operator()(dof_id_type dof_number,
                          const DofConstraintRow & constraint_row,

 * already constrains is skipped, as AddPrimalConstraint skips it, so that the
 * values collected are those the constraint path would have added.
 */
class CollectDirichletValues : public AddConstraint
{
private:
  const DofMap            & dof_map;
  DofConstraintValueMap   & values;

public:
  CollectDirichletValues(const DofMap & dof_map_in,
                         DofConstraintValueMap & values_in)
    : dof_map(dof_map_in), values(values_in) {}

  virtual void operator()(dof_id_type dof_number,
                          const DofConstraintRow & /*constraint_row*/,
                          const Number constraint_rhs) const
  {
    if (!dof_map.is_constrained_dof(dof_number))
      values[dof_number] = constraint_rhs;
  }
};

  } // apply_dirichlet_impl

public:
  ConstrainDirichlet (const DofMap & dof_map_in,
                      const MeshBase & mesh_in,
                      const Real time_in,
                      const DirichletBoundaries & dirichlets_in,

#ifdef LIBMESH_ENABLE_DIRICHLET
void DofMap::compute_dirichlet_values(const DirichletBoundaries & dirichlets,
                                      const MeshBase & mesh,
                                      const Real time,
                                      DofConstraintValueMap & values) const
{
  parallel_object_only();

  values.clear();

  if (dirichlets.empty())
    return;

  if (_verify_dirichlet_bc_consistency)
    for (const auto & dirichlet : dirichlets)
      this->check_dirichlet_bcid_consistency(mesh, *dirichlet);

  // Processors only project their local elements, as the constraint path does
  ConstElemRange range (mesh.local_elements_begin(),
                        mesh.local_elements_end());

  Threads::parallel_for
    (range, ConstrainDirichlet(*this, mesh, time, dirichlets,
                               CollectDirichletValues(*this, values)));
}
#endif // LIBMESH_ENABLE_DIRICHLET

      // vertices and scaled to unit H1 seminorm over the interval
    default:
      for (unsigned int n=1; n <= i; ++n)
        returnval *= xi;

      returnval = (returnval - ((i % 2) ? xi : 1.)) *
        fe_hierarchic_bubble_scaling(i);
      break;
    }

      // contributes, under the same scaling their shape functions carry
    default:
      for (unsigned int n=1; n != i; ++n)
        returnval *= xi;

      returnval = (Real(i) * returnval - ((i % 2) ? 1. : 0.)) *
        fe_hierarchic_bubble_scaling(i);
      break;
    }

      // Both parities differentiate twice to p (p-1) xi^(p-2), the linear term of an odd bubble
      // dropping out, under the same scaling their shape functions carry
    default:
      for (unsigned int n=2; n != i; ++n)
        returnval *= xi;

      returnval = Real(i) * (Real(i) - 1.) * returnval *
        fe_hierarchic_bubble_scaling(i);
      break;
    }

    {
      // The limit of the general expression below, in which only the bubble's leading term survives
      // and so carries the same normalization the one-dimensional bubble does
      return std::pow(edgenumerator, basisorder) *
        fe_hierarchic_bubble_scaling(basisorder);
    }
  // Experimentally, as c -> 0, n propto c, I'm still seeing good
  // behavior from the default implementation below:

      // Rotating the least node of the side to the origin leaves the
      // reflection about the diagonal through it, which is settled by
      // which of that node's two neighbors is the lesser.
      const bool flip = (side.point((min_side_node+3)%4) <
                         side.point((min_side_node+1)%4));

      switch (min_side_node) {
      case 0:

      case 1:
        sidep(0) = -sidep(0);
        if (!flip)
          std::swap(sidep(0), sidep(1));
        break;
      case 2:
        sidep(0) = -sidep(0);

      case 3:
        sidep(1) = -sidep(1);
        if (!flip)
          std::swap(sidep(0), sidep(1));
        break;
      default:
        libmesh_error();

      i1 = 0;
      i2 = 0;
      if (elem->positive_edge_orientation(0))
        xi = -xi_saved;
    }
  // Edge 1
  else if (i < 8 + 2*e)

      i1 = 1;
      i2 = 0;
      if (!elem->positive_edge_orientation(2))
        xi = -xi_saved;
    }
  // Edge 3
  else if (i < 8 + 4*e)

      i1 = 0;
      i2 = 1;
      if (elem->positive_edge_orientation(8))
        xi = -xi_saved;
    }
  // Edge 9
  else if (i < 8 + 10*e)

      i1 = 1;
      i2 = 1;
      if (!elem->positive_edge_orientation(10))
        xi = -xi_saved;
    }
  // Edge 11
  else if (i < 8 + 12*e)

        if (elem->positive_face_orientation(0))
          {
            // Case 1
            xi  = xi_saved;
            eta = eta_saved;
          }
        else
          {

        if (elem->positive_face_orientation(0))
          {
            // Case 3
            xi  = -eta_saved;
            eta = xi_saved;
          }
        else
          {

            eta = -eta_saved;
          }

      else if (elem->point(2) == min_point)
        if (elem->positive_face_orientation(0))
          {
            // Case 5
            xi  = -xi_saved;
            eta = -eta_saved;
          }
        else
          {
            // Case 6
            xi  = -eta_saved;
            eta = -xi_saved;
          }

      else if (elem->point(1) == min_point)
        {
          if (elem->positive_face_orientation(0))
            {
              // Case 7
              xi  = eta_saved;
              eta = -xi_saved;
            }
          else
            {
              // Case 8
              xi  = -xi_saved;
              eta = eta_saved;
            }
        }
    }

        if (!elem->positive_face_orientation(1))
          {
            // Case 1
            xi   = xi_saved;
            zeta = zeta_saved;
          }
        else
          {

          }

      else if (elem->point(1) == min_point)
        if (!elem->positive_face_orientation(1))
          {
            // Case 3
            xi   = zeta_saved;
            zeta = -xi_saved;
          }
        else
          {
            // Case 4
            xi   = -xi_saved;
            zeta = zeta_saved;
          }

      else if (elem->point(5) == min_point)
        if (!elem->positive_face_orientation(1))
          {
            // Case 5
            xi   = -xi_saved;
            zeta = -zeta_saved;
          }
        else
          {
            // Case 6
            xi   = -zeta_saved;
            zeta = -xi_saved;
          }

      else if (elem->point(4) == min_point)

          if (!elem->positive_face_orientation(1))
            {
              // Case 7
              xi   = -zeta_saved;
              zeta = xi_saved;
            }
          else
            {
              // Case 8
              xi   = xi_saved;
              zeta = -zeta_saved;
            }
        }
    }

        if (!elem->positive_face_orientation(2))
          {
            // Case 1
            eta  = eta_saved;
            zeta = zeta_saved;
          }
        else
          {

          }

      else if (elem->point(2) == min_point)
        if (!elem->positive_face_orientation(2))
          {
            // Case 3
            eta  = zeta_saved;
            zeta = -eta_saved;
          }
        else
          {
            // Case 4
            eta  = -eta_saved;
            zeta = zeta_saved;
          }

      else if (elem->point(6) == min_point)
        if (!elem->positive_face_orientation(2))
          {
            // Case 5
            eta  = -eta_saved;
            zeta = -zeta_saved;
          }
        else
          {

            zeta = -eta_saved;
          }

      else if (elem->point(5) == min_point)
        {
          if (!elem->positive_face_orientation(2))
            {
              // Case 7
              eta  = -zeta_saved;
              zeta = eta_saved;
            }
          else
            {
              // Case 8
              eta   = eta_saved;
              zeta = -zeta_saved;
            }
        }
    }

        if (elem->positive_face_orientation(3))
          {
            // Case 1
            xi   = xi_saved;
            zeta = zeta_saved;
          }
        else
          {

        else
          {
            // Case 4
            xi   = xi_saved;
            zeta = -zeta_saved;
          }

      else if (elem->point(6) == min_point)
        if (elem->positive_face_orientation(3))
          {
            // Case 5
            xi   = -xi_saved;
            zeta = -zeta_saved;
          }
        else
          {
            // Case 6
            xi   = -zeta_saved;
            zeta = -xi_saved;
          }

      else if (elem->point(2) == min_point)
        {
          if (elem->positive_face_orientation(3))
            {
              // Case 7
              xi   = zeta_saved;
              zeta = -xi_saved;
            }
          else
            {
              // Case 8
              xi   = -xi_saved;
              zeta = zeta_saved;
            }
        }
    }

        if (elem->positive_face_orientation(4))
          {
            // Case 1
            eta  = eta_saved;
            zeta = zeta_saved;
          }
        else
          {

          }

      else if (elem->point(4) == min_point)
        if (elem->positive_face_orientation(4))
          {
            // Case 3
            eta  = -zeta_saved;
            zeta = eta_saved;
          }
        else
          {
            // Case 4
            eta  = eta_saved;
            zeta = -zeta_saved;
          }

      else if (elem->point(7) == min_point)
        if (elem->positive_face_orientation(4))
          {
            // Case 5
            eta  = -eta_saved;
            zeta = -zeta_saved;
          }
        else
          {

            zeta = -eta_saved;
          }

      else if (elem->point(3) == min_point)
        {
          if (elem->positive_face_orientation(4))
            {
              // Case 7
              eta   = zeta_saved;
              zeta = -eta_saved;
            }
          else
            {
              // Case 8
              eta  = -eta_saved;
              zeta = zeta_saved;
            }
        }
    }

        if (!elem->positive_face_orientation(5))
          {
            // Case 1
            xi  = xi_saved;
            eta = eta_saved;
          }
        else
          {

          }

      else if (elem->point(5) == min_point)
        if (!elem->positive_face_orientation(5))
          {
            // Case 3
            xi  = eta_saved;
            eta = -xi_saved;
          }
        else
          {
            // Case 4
            xi  = -xi_saved;
            eta = eta_saved;
          }

      else if (elem->point(6) == min_point)
        if (!elem->positive_face_orientation(5))
          {
            // Case 5
            xi  = -xi_saved;
            eta = -eta_saved;
          }
        else
          {
            // Case 6
            xi  = -eta_saved;
            eta = -xi_saved;
          }

      else if (elem->point(7) == min_point)

          if (!elem->positive_face_orientation(5))
            {
              // Case 7
              xi  = -eta_saved;
              eta = xi_saved;
            }
          else
            {
              // Case 8
              xi  = xi_saved;
              eta = -eta_saved;
            }
        }
    }

      i1 = cube_number_row[basisnum] + 2;
      i2 = cube_number_page[basisnum] + 2;
    }
}


// Reorder the barycentric coordinates of a triangular face of a prism, whose vertices begin at
// \p first_vertex, so that they follow the order of the face's vertices. The interior basis of a
// triangle is not symmetric in its barycentric coordinates, so a face shared between two elements
// needs them ordered the same way from both sides.
void orient_triangle_coords(const Elem & elem,
                            const unsigned int first_vertex,
                            const Point & xi_eta_saved,
                            Point & xi_eta)
{
  unsigned int face_vertex[3] = {first_vertex, first_vertex+1, first_vertex+2};
  orient_triangle(elem, face_vertex);

  const Real barycentric[3] = {1 - xi_eta_saved(0) - xi_eta_saved(1),
                               xi_eta_saved(0),
                               xi_eta_saved(1)};

  xi_eta(0) = barycentric[face_vertex[1] - first_vertex];
  xi_eta(1) = barycentric[face_vertex[2] - first_vertex];
}

      if (elem->point(0) == min_point)
        {
          if (!elem->positive_face_orientation(1))
            {
              // Case 1: no flips needed
              i01 = s0+1; // edge to triangle side 0 numbering

          else
            {
              // Case 2: flip about 0-4 diagonal
              i01 = s1+1;
              i2 = s0;
            }
        }
      else if (elem->point(3) == min_point)

              // Case 3: 0->3->4->1->0 rotation
              i01 = s1+1;
              i2 = s0;
              zeta = -zeta_saved;
            }
          else
            {

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

              // Case 8: flip about 1-3 diagonal
              i01 = s1+1;
              i2 = s0;
              xi_eta(0) = (1-xe_fraction)*xe_scale;
              zeta = -zeta_saved;
            }
        }
    }

      if (elem->point(1) == min_point)
        {
          if (!elem->positive_face_orientation(2))
            {
              // Case 1: no flips needed
              i01 = s0+1+e; // edge to triangle side 1 numbering
              i2 = s1;
            }
          else
            {
              // Case 2: flip about 1-5 diagonal
              i01 = s1+1+e;
              i2 = s0;
            }
        }
      else if (elem->point(4) == min_point)

              // Case 3: 1->4->5->2->1 rotation
              i01 = s1+1+e;
              i2 = s0;
              zeta = -zeta_saved;
            }
          else
            {

              zeta = -zeta_saved;
            }
        }
      else if (elem->point(2) == min_point)
        {
          if (!elem->positive_face_orientation(2))
            {
              // Case 5: 1->2->5->4->1 rotation
              i01 = s1+1+e;
              i2 = s0;
              const Real xe = xe_fraction;
              xi_eta(1) = xe*xe_scale;
              xi_eta(0) = xe_scale - xi_eta(1);
            }
          else
            {
              // Case 6: flip about 7-13 midline
              i01 = s0+1+e;
              i2 = s1;
              const Real xe = xe_fraction;
              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(2))
            {

              i01 = s0+1+e;
              i2 = s1;
              zeta = -zeta_saved;
              const Real xe = xe_fraction;
              xi_eta(1) = xe*xe_scale;
              xi_eta(0) = xe_scale - xi_eta(1);
            }

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

      if (elem->point(2) == min_point)
        {
          if (!elem->positive_face_orientation(3))
            {
              // Case 1: no flips needed
              i01 = s0+1+2*e; // edge to triangle side 2 numbering

          else
            {
              // Case 2: flip about 2-3 diagonal
              i01 = s1+1+2*e;
              i2 = s0;
            }
        }
      else if (elem->point(5) == min_point)

              // Case 3: 2->5->3->0->2 rotation
              i01 = s1+1+2*e;
              i2 = s0;
              zeta = -zeta_saved;
            }
          else
            {

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

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

    {
      i01 = i - 3 - 6*e - 3*e*e;
      i2 = 0;
      orient_triangle_coords(*elem, 0, xi_eta_saved, xi_eta);
    }
  // Face 4
  else if (i < 6 + 9*e + 3*e*e + e*(e-1))
    {
      i01 = i - 3 - 6*e - 3*e*e - e*(e-1)/2;
      i2 = 1;
      orient_triangle_coords(*elem, 3, xi_eta_saved, xi_eta);
    }
  // Internal DoFs
  else

            if (i01 > 2 && i01 < 3u*totalorder)
              {
                // %(p-1) to find the edge number, %2 for even vs odd
                const bool odd_basis = ((i01-3)%(totalorder-1))%2;
                if (odd_basis)
                  {
                    const int tri_edge = (i01-3)/(totalorder-1);

              {
                // The limit of the general expression below, in which only the bubble's leading term
                // survives and so carries the same normalization the one-dimensional bubble does
                return std::pow(edgenumerator, basisorder) *
                  fe_hierarchic_bubble_scaling(basisorder);
              }

            const Real edgeval = edgenumerator / crossval;

bool
Elem::positive_face_orientation(const unsigned int i) const
{
  return this->face_orientation(i) % 2;
}


unsigned int
Elem::edge_orientation(const unsigned int i) const
{
  return this->positive_edge_orientation(i);
}


unsigned int
Elem::face_orientation(const unsigned int i) const
{
  libmesh_assert_less (i, this->n_faces());

  const unsigned int v = std::distance(nodes.begin(),
                         std::min_element(nodes.begin(), nodes.begin() + N, cmp));

  return 2 * v + cmp(nodes[(v - 1 + N) % N], nodes[(v + 1) % N]);
}

bool

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

    {
      ray_target = inside - Point(1);
      intersection_distances =
        this->find_ray_intersections(inside, ray_target);
    }

  // I'd make this an assert, but I'm not 100% confident we can't

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

                                const numeric_index_type,
                                const numeric_index_type blocksize)
{
  init_shell_mat(*this, m, n, m_l, n_l, blocksize);
  this->set_context();
}

template <typename T>
void

template <typename T>
void
PetscMatrixShellMatrix<T>::zero()
{
  // A shell matrix computes its action and stores no entries, so there is nothing to clear. This is
  // reachable through System::init_matrices(), which zeroes every matrix it initializes.
}

template <typename T>
std::unique_ptr<SparseMatrix<T>>
PetscMatrixShellMatrix<T>::zero_clone() const
{
  libmesh_error();
}

template <typename T>
std::unique_ptr<SparseMatrix<T>>
PetscMatrixShellMatrix<T>::clone() const
{
  libmesh_not_implemented();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::set(const numeric_index_type, const numeric_index_type, const T)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const numeric_index_type, const numeric_index_type, const T)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
                                      const std::vector<numeric_index_type> &,
                                      const std::vector<numeric_index_type> &)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add_matrix(const DenseMatrix<T> &,
                                      const std::vector<numeric_index_type> &)
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::add(const T, const SparseMatrix<T> &)
{
  libmesh_error();
}

template <typename T>
T
PetscMatrixShellMatrix<T>::operator()(const numeric_index_type, const numeric_index_type) const
{
  libmesh_error();
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::l1_norm() const
{
  libmesh_error();
}

template <typename T>
Real
PetscMatrixShellMatrix<T>::linfty_norm() const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::print_personal(std::ostream &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_diagonal(NumericVector<T> &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_transpose(SparseMatrix<T> &) const
{
  libmesh_error();
}

template <typename T>
void
PetscMatrixShellMatrix<T>::get_row(numeric_index_type,
                                   std::vector<numeric_index_type> &,
                                   std::vector<T> &) const
{
  libmesh_error();
}

template class LIBMESH_EXPORT PetscMatrixShellMatrix<Number>;

        Jac = &mffd_jac;
        // mffd_jac is function-local, so don't attach a context to jac here -- it would
        // dangle once mffd_jac is destroyed at the end of this call.
        mffd_jac.assign(jac, /*set_context=*/false);
      }

    // We already computed the Jacobian during the residual evaluation

  _default_monitor(true),
  _snesmf_reuse_base(true),
  _computing_base_vector(true),
  _setup_reuse(false)
{
}

                                const double        tol,     // Stopping tolerance
                                const unsigned int  m_its)
{
  return this->solve(pre_in, pre_in, x_in, r_in, tol, m_its);
}

template <typename T>
std::pair<unsigned int, Real>
PetscNonlinearSolver<T>::solve (SparseMatrix<T> &  jac_in,  // Jacobian operator matrix (Amat)
                                SparseMatrix<T> &  pre_in,  // Preconditioning matrix (Pmat)
                                NumericVector<T> & x_in,    // Solution vector
                                NumericVector<T> & r_in,    // Residual vector

  this->init ();

  // Make sure the data passed in are really of Petsc types
  PetscMatrixBase<T> * jac = cast_ptr<PetscMatrixBase<T> *>(&jac_in);
  PetscMatrixBase<T> * pre = cast_ptr<PetscMatrixBase<T> *>(&pre_in);
  PetscVector<T> * x   = cast_ptr<PetscVector<T> *>(&x_in);
  PetscVector<T> * r   = cast_ptr<PetscVector<T> *>(&r_in);

  // Only set the jacobian function if we've been provided with something to call.
  // This allows a user to set their own jacobian function if they want to
  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
    LibmeshPetscCall(SNESSetJacobian (_snes, jac->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));

  // Have the Krylov subspace method use our good initial guess rather than 0
  KSP ksp;

  nonlinear_solver          (NonlinearSolver<Number>::build(*this)),
  diff_solver               (),
  _n_nonlinear_iterations   (0),
  _final_nonlinear_residual (1.e20),
  _operator_matrix          (nullptr)
{
  // Set default parameters
  // These were chosen to match the Petsc defaults

    es.parameters.get<unsigned int>("nonlinear solver maximum function evaluations");

  const double abs_resid_tol = parameters.have_parameter<Real>("nonlinear solver absolute residual tolerance") ?
    double(parameters.get<Real>("nonlinear solver absolute residual tolerance")) :
    double(es.parameters.get<Real>("nonlinear solver absolute residual tolerance"));

  const double rel_resid_tol = parameters.have_parameter<Real>("nonlinear solver relative residual tolerance") ?

  else
    {
      if (this->prefix_with_name())
        nonlinear_solver->init(this->prefix().c_str());
      else
        nonlinear_solver->init();

      // If a distinct Jacobian operator matrix has been registered (see
      // set_operator_matrix()), use it as Amat while *matrix remains the preconditioning matrix
      // (Pmat); otherwise use the ordinary single-matrix solve.
      if (_operator_matrix)
        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
          nonlinear_solver->solve (*_operator_matrix, *matrix, *solution, *rhs,
                                   nonlinear_solver->relative_residual_tolerance,
                                   nonlinear_solver->max_linear_iterations);
      else
        std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
          nonlinear_solver->solve (*matrix, *solution, *rhs,
                                   nonlinear_solver->relative_residual_tolerance,
                                   nonlinear_solver->max_linear_iterations);
    }

  // Update the system after the solve