From 29d93763c3adf5247340568d9f8582dd387ad7b2 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 10:56:27 +0100
Subject: [PATCH 1/7] Use userspace solver

---
 cpp/demo/hyperelasticity/CMakeLists.txt     |   2 +-
 cpp/demo/hyperelasticity/NewtonSolver.cpp   | 199 ++++++++++++++++++++
 cpp/demo/hyperelasticity/NewtonSolver.h     | 190 +++++++++++++++++++
 cpp/demo/hyperelasticity/hyperelasticity.py |   7 +-
 cpp/demo/hyperelasticity/main.cpp           |   6 +-
 5 files changed, 398 insertions(+), 6 deletions(-)
 create mode 100644 cpp/demo/hyperelasticity/NewtonSolver.cpp
 create mode 100644 cpp/demo/hyperelasticity/NewtonSolver.h

diff --git a/cpp/demo/hyperelasticity/CMakeLists.txt b/cpp/demo/hyperelasticity/CMakeLists.txt
index 007b04e24c5..77273759b1f 100644
--- a/cpp/demo/hyperelasticity/CMakeLists.txt
+++ b/cpp/demo/hyperelasticity/CMakeLists.txt
@@ -35,7 +35,7 @@ else()
 
   set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-  add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
+  add_executable(${PROJECT_NAME} main.cpp NewtonSolver.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
   target_link_libraries(${PROJECT_NAME} dolfinx)
 
   # Set C++20 standard
diff --git a/cpp/demo/hyperelasticity/NewtonSolver.cpp b/cpp/demo/hyperelasticity/NewtonSolver.cpp
new file mode 100644
index 00000000000..186687033c0
--- /dev/null
+++ b/cpp/demo/hyperelasticity/NewtonSolver.cpp
@@ -0,0 +1,199 @@
+// Copyright (C) 2005-2021 Garth N. Wells
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+// #ifdef HAS_PETSC
+
+#include "NewtonSolver.h"
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/common/log.h>
+#include <dolfinx/la/petsc.h>
+#include <string>
+#include <utility>
+
+using namespace dolfinx;
+
+/// Convergence test
+/// @param solver The Newton solver
+/// @param r The residual vector
+/// @return The pair `(residual norm, converged)`, where `converged` is
+/// and true` if convergence achieved
+std::pair<double, bool> NewtonSolver::converged(const Vec r)
+{
+  PetscReal residual = 0;
+  VecNorm(r, NORM_2, &residual);
+
+  // Relative residual
+  const double relative_residual = residual / residual0();
+
+  // Output iteration number and residual
+  if (dolfinx::MPI::rank(comm()) == 0)
+  {
+    spdlog::info("Newton iteration {}"
+                 ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
+                 _iteration, residual, atol, relative_residual, rtol);
+  }
+
+  // Return true if convergence criterion is met
+  if (relative_residual < rtol or residual < atol)
+    return {residual, true};
+  else
+    return {residual, false};
+}
+
+//-----------------------------------------------------------------------------
+NewtonSolver::NewtonSolver(MPI_Comm comm)
+    : _krylov_iterations(0), _iteration(0), _residual(0), _residual0(0),
+      _solver(comm), _dx(nullptr), _comm(comm)
+{
+  // Create linear solver if not already created. Default to LU.
+  _solver.set_options_prefix("nls_solve_");
+  la::petsc::options::set("nls_solve_ksp_type", "preonly");
+  la::petsc::options::set("nls_solve_pc_type", "lu");
+  _solver.set_from_options();
+}
+//-----------------------------------------------------------------------------
+NewtonSolver::~NewtonSolver()
+{
+  if (_b)
+    VecDestroy(&_b);
+  if (_dx)
+    VecDestroy(&_dx);
+  if (_matJ)
+    MatDestroy(&_matJ);
+  if (_matP)
+    MatDestroy(&_matP);
+}
+//-----------------------------------------------------------------------------
+void NewtonSolver::setF(std::function<void(const Vec, Vec)> F,
+                                    Vec b)
+{
+  _fnF = std::move(F);
+  _b = b;
+  PetscObjectReference((PetscObject)_b);
+}
+//-----------------------------------------------------------------------------
+void NewtonSolver::setJ(std::function<void(const Vec, Mat)> J,
+                                    Mat Jmat)
+{
+  _fnJ = std::move(J);
+  _matJ = Jmat;
+  PetscObjectReference((PetscObject)_matJ);
+}
+//-----------------------------------------------------------------------------
+const la::petsc::KrylovSolver&
+NewtonSolver::get_krylov_solver() const
+{
+  return _solver;
+}
+//-----------------------------------------------------------------------------
+la::petsc::KrylovSolver& NewtonSolver::get_krylov_solver()
+{
+  return _solver;
+}
+//-----------------------------------------------------------------------------
+void NewtonSolver::set_form(std::function<void(Vec)> form)
+{
+  _system = std::move(form);
+}
+//-----------------------------------------------------------------------------
+std::pair<int, bool> NewtonSolver::solve(Vec x)
+{
+  // Reset iteration counts
+  _iteration = 0;
+  int krylov_iterations = 0;
+  _residual = -1;
+  _residual0 = 0;
+
+  if (!_fnF)
+  {
+    throw std::runtime_error("Function for computing residual vector has not "
+                             "been provided to the NewtonSolver.");
+  }
+
+  if (!_fnJ)
+  {
+    throw std::runtime_error("Function for computing Jacobian has not "
+                             "been provided to the NewtonSolver.");
+  }
+
+  if (_system)
+    _system(x);
+  assert(_b);
+  _fnF(x, _b);
+
+  // Check convergence
+  bool newton_converged = false;
+  if (convergence_criterion == "residual")
+    std::tie(_residual, newton_converged) = converged(_b);
+
+  _solver.set_operators(_matJ, _matJ);
+
+  if (!_dx)
+    MatCreateVecs(_matJ, &_dx, nullptr);
+
+  // Start iterations
+  while (!newton_converged and _iteration < max_it)
+  {
+    // Compute Jacobian
+    assert(_matJ);
+    _fnJ(x, _matJ);
+
+    // Perform linear solve and update total number of Krylov iterations
+    krylov_iterations += _solver.solve(_dx, _b);
+
+    // Update solution
+    VecAXPY(x, -relaxation_parameter, _dx);
+
+    // Increment iteration count
+    ++_iteration;
+
+    // Compute F
+    if (_system)
+      _system(x);
+    _fnF(x, _b);
+    // Initialize _residual0
+    if (_iteration == 1)
+    {
+      PetscReal _r = 0;
+      VecNorm(_dx, NORM_2, &_r);
+      _residual0 = _r;
+    }
+
+    // Test for convergence
+    std::tie(_residual, newton_converged) = converged(_b);
+  }
+
+  if (newton_converged)
+  {
+    if (dolfinx::MPI::rank(_comm.comm()) == 0)
+    {
+      spdlog::info("Newton solver finished in {} iterations and {} linear "
+                   "solver iterations.",
+                   _iteration, krylov_iterations);
+    }
+  }
+  else
+  {
+    throw std::runtime_error("Newton solver did not converge.");
+  }
+
+  return {_iteration, newton_converged};
+}
+//-----------------------------------------------------------------------------
+int NewtonSolver::krylov_iterations() const
+{
+  return _krylov_iterations;
+}
+//-----------------------------------------------------------------------------
+// int NewtonSolver::iteration() const { return _iteration; }
+//-----------------------------------------------------------------------------
+// double NewtonSolver::residual() const { return _residual; }
+//-----------------------------------------------------------------------------
+double NewtonSolver::residual0() const { return _residual0; }
+//-----------------------------------------------------------------------------
+MPI_Comm NewtonSolver::comm() const { return _comm.comm(); }
+//-----------------------------------------------------------------------------
+// #endif
diff --git a/cpp/demo/hyperelasticity/NewtonSolver.h b/cpp/demo/hyperelasticity/NewtonSolver.h
new file mode 100644
index 00000000000..70907c2537e
--- /dev/null
+++ b/cpp/demo/hyperelasticity/NewtonSolver.h
@@ -0,0 +1,190 @@
+// Copyright (C) 2005-2021 Garth N. Wells
+//
+// This file is part of DOLFINx (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#pragma once
+
+// #ifdef HAS_PETSC
+
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/la/petsc.h>
+#include <functional>
+#include <memory>
+#include <petscmat.h>
+#include <petscvec.h>
+#include <utility>
+
+
+/// @brief A Newton solver for nonlinear systems of equations of the
+/// form \f$F(x) = 0\f$.
+///
+/// It solves \f[ \left. \frac{dF}{dx} \right|_{x} \Delta x = F(x) \f]
+/// with default update \f$x \leftarrow x - \Delta x\f$.
+///
+/// It relies on PETSc for linear algebra backends.
+class NewtonSolver
+{
+public:
+  /// @brief Create nonlinear solver
+  /// @param[in] comm MPI communicator for the solver
+  explicit NewtonSolver(MPI_Comm comm);
+
+  /// Move constructor
+  NewtonSolver(NewtonSolver&& solver) = default;
+
+  // Copy constructor (deleted)
+  NewtonSolver(const NewtonSolver& solver) = delete;
+
+  /// Move assignment constructor
+  NewtonSolver& operator=(NewtonSolver&& solver) = default;
+
+  // Assignment operator (deleted)
+  NewtonSolver& operator=(const NewtonSolver& solver) = delete;
+
+  /// Destructor
+  ~NewtonSolver();
+
+  std::pair<double, bool> converged(const Vec r);
+
+  /// @brief Set the function for computing the residual \f$F(x) = 0\f$
+  /// and the vector to assemble the residual into.
+  /// @param[in] F Function to compute/assemble the residual vector `b`.
+  /// The first argument to the function is the solution vector `x` and
+  /// the second is the vector `b` to assemble into.
+  /// @param[in] b Vector to assemble to residual into.
+  void setF(std::function<void(const Vec, Vec)> F, Vec b);
+
+  /// @brief Set the function for computing the Jacobian \f$J:=dF/dx\f$
+  /// and the matrix to assemble the Jacobian into.
+  /// @param[in] J Function to compute the Jacobian matrix `Jmat`. The
+  /// first argument to the function is the solution vector `x` and the
+  /// second is the matrix to assemble into.
+  /// @param[in] Jmat Matrix to assemble the Jacobian into.
+  void setJ(std::function<void(const Vec, Mat)> J, Mat Jmat);
+
+  /// @brief Get the internal Krylov solver used to solve for the Newton
+  /// updates (const version).
+  ///
+  /// The Krylov solver prefix is `nls_solve_`.
+  ///
+  /// @return The Krylov solver
+  const dolfinx::la::petsc::KrylovSolver& get_krylov_solver() const;
+
+  /// @brief Get the internal Krylov solver used to solve for the Newton
+  /// updates (non-const version).
+  ///
+  /// The Krylov solver prefix is `nls_solve_`.
+  ///
+  /// @return The Krylov solver
+  dolfinx::la::petsc::KrylovSolver& get_krylov_solver();
+
+  /// @brief Set the function that is called before the residual or
+  /// Jacobian are computed. It is commonly used to update ghost values.
+  ///
+  /// @param[in] form Function to call. It takes the (latest) solution
+  /// vector `x` as an argument.
+  void set_form(std::function<void(Vec)> form);
+
+  /// @brief Optional set function that is called after each inner solve for
+  /// the Newton increment to update the solution.
+  ///
+  /// The function `update` takes `this`, the Newton increment `dx`, and
+  /// the vector `x` from the start of the Newton iteration.
+  ///
+  /// By default, the update is x <- x - dx
+  ///
+  /// @param[in] update The function that updates the solution
+  void set_update(
+      std::function<void(const NewtonSolver& solver, const Vec, Vec)> update);
+
+  /// @brief Solve the nonlinear problem.
+  ///
+  /// @param[in,out] x The solution vector. It should be set the initial
+  /// solution guess.
+  /// @return (number of Newton iterations, whether iteration converged)
+  std::pair<int, bool> solve(Vec x);
+
+  /// @brief Number of Krylov iterations elapsed since solve started.
+  /// @return Number of Krylov iterations.
+  int krylov_iterations() const;
+
+  /// @brief Get initial residual.
+  /// @return Initial residual.
+  double residual0() const;
+
+  /// @brief Get MPI communicator.
+  MPI_Comm comm() const;
+
+  /// @brief Maximum number of iterations.
+  int max_it = 50;
+
+  /// @brief Relative convergence tolerance.
+  double rtol = 1e-9;
+
+  /// @brief Absolute convergence tolerance.
+  double atol = 1e-10;
+
+  /// @todo change to string to enum.
+  /// @brief Convergence criterion.
+  std::string convergence_criterion = "residual";
+
+  /// @brief Relaxation parameter.
+  double relaxation_parameter = 1.0;
+
+private:
+  // Function for computing the residual vector. The first argument is
+  // the latest solution vector x and the second argument is the
+  // residual vector.
+  std::function<void(const Vec x, Vec b)> _fnF;
+
+  // Function for computing the Jacobian matrix operator. The first
+  // argument is the latest solution vector x and the second argument is
+  // the matrix operator.
+  std::function<void(const Vec x, Mat J)> _fnJ;
+
+  // Function for computing the preconditioner matrix operator. The
+  // first argument is the latest solution vector x and the second
+  // argument is the matrix operator.
+  std::function<void(const Vec x, Mat P)> _fnP;
+
+  // Function called before the residual and Jacobian function at each
+  // iteration.
+  std::function<void(const Vec x)> _system;
+
+  // Residual vector
+  Vec _b = nullptr;
+
+  // Jacobian matrix and preconditioner matrix
+  Mat _matJ = nullptr, _matP = nullptr;
+
+  // Function to check for convergence
+  std::function<std::pair<double, bool>(const NewtonSolver& solver,
+                                        const Vec r)>
+      _converged;
+
+  // Function to update the solution after inner solve for the Newton increment
+  std::function<void(const NewtonSolver& solver, const Vec dx, Vec x)>
+      _update_solution;
+
+  // Accumulated number of Krylov iterations since solve began
+  int _krylov_iterations;
+
+  // Number of iterations
+  int _iteration;
+
+  // Most recent residual and initial residual
+  double _residual, _residual0;
+
+  // Linear solver
+  dolfinx::la::petsc::KrylovSolver _solver;
+
+  // Solution vector
+  Vec _dx = nullptr;
+
+  // MPI communicator
+  dolfinx::MPI::Comm _comm;
+};
+
+// #endif
diff --git a/cpp/demo/hyperelasticity/hyperelasticity.py b/cpp/demo/hyperelasticity/hyperelasticity.py
index f361c4171ab..97d6d6ccd6e 100644
--- a/cpp/demo/hyperelasticity/hyperelasticity.py
+++ b/cpp/demo/hyperelasticity/hyperelasticity.py
@@ -10,6 +10,7 @@
     Constant,
     FunctionSpace,
     Identity,
+    Measure,
     Mesh,
     TestFunction,
     TrialFunction,
@@ -17,7 +18,6 @@
     det,
     diff,
     ds,
-    dx,
     grad,
     inner,
     ln,
@@ -27,8 +27,9 @@
 
 # Function spaces
 e = element("Lagrange", "tetrahedron", 1, shape=(3,))
+e2 = element("Lagrange", "tetrahedron", 2, shape=(3,))
 mesh = Mesh(e)
-V = FunctionSpace(mesh, e)
+V = FunctionSpace(mesh, e2)
 
 # Trial and test functions
 du = TrialFunction(V)  # Incremental displacement
@@ -74,6 +75,8 @@
 # Stored strain energy density (compressible neo-Hookean model)
 psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
 
+
+dx = Measure("dx", metadata={"quadrature_degree": 3})
 # Total potential energy
 Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds
 
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index df95ae8b51b..a4d35edcc61 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -14,6 +14,7 @@
 
 // ## C++ program
 
+#include "NewtonSolver.h"
 #include "hyperelasticity.h"
 #include <algorithm>
 #include <basix/finite-element.h>
@@ -28,7 +29,6 @@
 #include <dolfinx/la/petsc.h>
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/cell_types.h>
-#include <dolfinx/nls/NewtonSolver.h>
 #include <numbers>
 #include <petscmat.h>
 #include <petscsys.h>
@@ -162,7 +162,7 @@ int main(int argc, char* argv[])
         mesh::create_cell_partitioner(mesh::GhostMode::none)));
 
     auto element = basix::create_element<U>(
-        basix::element::family::P, basix::cell::type::tetrahedron, 1,
+        basix::element::family::P, basix::cell::type::tetrahedron, 2,
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
@@ -247,7 +247,7 @@ int main(int argc, char* argv[])
            fem::DirichletBC<T>(u_rotation, bdofs_right)};
 
     HyperElasticProblem problem(L, a, bcs);
-    nls::petsc::NewtonSolver newton_solver(mesh->comm());
+    NewtonSolver newton_solver(mesh->comm());
     newton_solver.setF(problem.F(), problem.vector());
     newton_solver.setJ(problem.J(), problem.matrix());
     newton_solver.set_form(problem.form());

From 484582674e7b9e1ea98911492ac9b0c975799a49 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 11:20:30 +0100
Subject: [PATCH 2/7] More trimming

---
 cpp/demo/hyperelasticity/NewtonSolver.cpp   | 107 +++++++-------------
 cpp/demo/hyperelasticity/NewtonSolver.h     |  40 +-------
 cpp/demo/hyperelasticity/hyperelasticity.py |   6 +-
 cpp/demo/hyperelasticity/main.cpp           |   2 +-
 4 files changed, 44 insertions(+), 111 deletions(-)

diff --git a/cpp/demo/hyperelasticity/NewtonSolver.cpp b/cpp/demo/hyperelasticity/NewtonSolver.cpp
index 186687033c0..3e76388d285 100644
--- a/cpp/demo/hyperelasticity/NewtonSolver.cpp
+++ b/cpp/demo/hyperelasticity/NewtonSolver.cpp
@@ -15,38 +15,9 @@
 
 using namespace dolfinx;
 
-/// Convergence test
-/// @param solver The Newton solver
-/// @param r The residual vector
-/// @return The pair `(residual norm, converged)`, where `converged` is
-/// and true` if convergence achieved
-std::pair<double, bool> NewtonSolver::converged(const Vec r)
-{
-  PetscReal residual = 0;
-  VecNorm(r, NORM_2, &residual);
-
-  // Relative residual
-  const double relative_residual = residual / residual0();
-
-  // Output iteration number and residual
-  if (dolfinx::MPI::rank(comm()) == 0)
-  {
-    spdlog::info("Newton iteration {}"
-                 ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
-                 _iteration, residual, atol, relative_residual, rtol);
-  }
-
-  // Return true if convergence criterion is met
-  if (relative_residual < rtol or residual < atol)
-    return {residual, true};
-  else
-    return {residual, false};
-}
-
 //-----------------------------------------------------------------------------
 NewtonSolver::NewtonSolver(MPI_Comm comm)
-    : _krylov_iterations(0), _iteration(0), _residual(0), _residual0(0),
-      _solver(comm), _dx(nullptr), _comm(comm)
+    : _solver(comm), _dx(nullptr), _comm(comm)
 {
   // Create linear solver if not already created. Default to LU.
   _solver.set_options_prefix("nls_solve_");
@@ -83,17 +54,6 @@ void NewtonSolver::setJ(std::function<void(const Vec, Mat)> J,
   PetscObjectReference((PetscObject)_matJ);
 }
 //-----------------------------------------------------------------------------
-const la::petsc::KrylovSolver&
-NewtonSolver::get_krylov_solver() const
-{
-  return _solver;
-}
-//-----------------------------------------------------------------------------
-la::petsc::KrylovSolver& NewtonSolver::get_krylov_solver()
-{
-  return _solver;
-}
-//-----------------------------------------------------------------------------
 void NewtonSolver::set_form(std::function<void(Vec)> form)
 {
   _system = std::move(form);
@@ -102,10 +62,34 @@ void NewtonSolver::set_form(std::function<void(Vec)> form)
 std::pair<int, bool> NewtonSolver::solve(Vec x)
 {
   // Reset iteration counts
-  _iteration = 0;
+  int iteration = 0;
   int krylov_iterations = 0;
-  _residual = -1;
-  _residual0 = 0;
+  double residual = -1;
+  double residual0 = 0;
+
+  auto converged = [&iteration, &residual, &residual0,
+                    this](const Vec r) -> std::pair<double, bool>
+  {
+    PetscReal residual = 0;
+    VecNorm(r, NORM_2, &residual);
+
+    // Relative residual
+    const double relative_residual = residual / residual0;
+
+    // Output iteration number and residual
+    if (dolfinx::MPI::rank(_comm.comm()) == 0)
+    {
+      spdlog::info("Newton iteration {}"
+                   ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
+                   iteration, residual, atol, relative_residual, rtol);
+    }
+
+    // Return true if convergence criterion is met
+    if (relative_residual < rtol or residual < atol)
+      return {residual, true};
+    else
+      return {residual, false};
+  };
 
   if (!_fnF)
   {
@@ -126,16 +110,14 @@ std::pair<int, bool> NewtonSolver::solve(Vec x)
 
   // Check convergence
   bool newton_converged = false;
-  if (convergence_criterion == "residual")
-    std::tie(_residual, newton_converged) = converged(_b);
+  std::tie(residual, newton_converged) = converged(_b);
 
   _solver.set_operators(_matJ, _matJ);
 
-  if (!_dx)
-    MatCreateVecs(_matJ, &_dx, nullptr);
+  MatCreateVecs(_matJ, &_dx, nullptr);
 
   // Start iterations
-  while (!newton_converged and _iteration < max_it)
+  while (!newton_converged and iteration < max_it)
   {
     // Compute Jacobian
     assert(_matJ);
@@ -148,22 +130,23 @@ std::pair<int, bool> NewtonSolver::solve(Vec x)
     VecAXPY(x, -relaxation_parameter, _dx);
 
     // Increment iteration count
-    ++_iteration;
+    ++iteration;
 
     // Compute F
     if (_system)
       _system(x);
     _fnF(x, _b);
+
     // Initialize _residual0
-    if (_iteration == 1)
+    if (iteration == 1)
     {
       PetscReal _r = 0;
       VecNorm(_dx, NORM_2, &_r);
-      _residual0 = _r;
+      residual0 = _r;
     }
 
     // Test for convergence
-    std::tie(_residual, newton_converged) = converged(_b);
+    std::tie(residual, newton_converged) = converged(_b);
   }
 
   if (newton_converged)
@@ -172,7 +155,7 @@ std::pair<int, bool> NewtonSolver::solve(Vec x)
     {
       spdlog::info("Newton solver finished in {} iterations and {} linear "
                    "solver iterations.",
-                   _iteration, krylov_iterations);
+                   iteration, krylov_iterations);
     }
   }
   else
@@ -180,20 +163,6 @@ std::pair<int, bool> NewtonSolver::solve(Vec x)
     throw std::runtime_error("Newton solver did not converge.");
   }
 
-  return {_iteration, newton_converged};
+  return {iteration, newton_converged};
 }
-//-----------------------------------------------------------------------------
-int NewtonSolver::krylov_iterations() const
-{
-  return _krylov_iterations;
-}
-//-----------------------------------------------------------------------------
-// int NewtonSolver::iteration() const { return _iteration; }
-//-----------------------------------------------------------------------------
-// double NewtonSolver::residual() const { return _residual; }
-//-----------------------------------------------------------------------------
-double NewtonSolver::residual0() const { return _residual0; }
-//-----------------------------------------------------------------------------
-MPI_Comm NewtonSolver::comm() const { return _comm.comm(); }
-//-----------------------------------------------------------------------------
 // #endif
diff --git a/cpp/demo/hyperelasticity/NewtonSolver.h b/cpp/demo/hyperelasticity/NewtonSolver.h
index 70907c2537e..82d6ff2811b 100644
--- a/cpp/demo/hyperelasticity/NewtonSolver.h
+++ b/cpp/demo/hyperelasticity/NewtonSolver.h
@@ -46,8 +46,6 @@ class NewtonSolver
   /// Destructor
   ~NewtonSolver();
 
-  std::pair<double, bool> converged(const Vec r);
-
   /// @brief Set the function for computing the residual \f$F(x) = 0\f$
   /// and the vector to assemble the residual into.
   /// @param[in] F Function to compute/assemble the residual vector `b`.
@@ -64,22 +62,6 @@ class NewtonSolver
   /// @param[in] Jmat Matrix to assemble the Jacobian into.
   void setJ(std::function<void(const Vec, Mat)> J, Mat Jmat);
 
-  /// @brief Get the internal Krylov solver used to solve for the Newton
-  /// updates (const version).
-  ///
-  /// The Krylov solver prefix is `nls_solve_`.
-  ///
-  /// @return The Krylov solver
-  const dolfinx::la::petsc::KrylovSolver& get_krylov_solver() const;
-
-  /// @brief Get the internal Krylov solver used to solve for the Newton
-  /// updates (non-const version).
-  ///
-  /// The Krylov solver prefix is `nls_solve_`.
-  ///
-  /// @return The Krylov solver
-  dolfinx::la::petsc::KrylovSolver& get_krylov_solver();
-
   /// @brief Set the function that is called before the residual or
   /// Jacobian are computed. It is commonly used to update ghost values.
   ///
@@ -106,17 +88,6 @@ class NewtonSolver
   /// @return (number of Newton iterations, whether iteration converged)
   std::pair<int, bool> solve(Vec x);
 
-  /// @brief Number of Krylov iterations elapsed since solve started.
-  /// @return Number of Krylov iterations.
-  int krylov_iterations() const;
-
-  /// @brief Get initial residual.
-  /// @return Initial residual.
-  double residual0() const;
-
-  /// @brief Get MPI communicator.
-  MPI_Comm comm() const;
-
   /// @brief Maximum number of iterations.
   int max_it = 50;
 
@@ -126,10 +97,6 @@ class NewtonSolver
   /// @brief Absolute convergence tolerance.
   double atol = 1e-10;
 
-  /// @todo change to string to enum.
-  /// @brief Convergence criterion.
-  std::string convergence_criterion = "residual";
-
   /// @brief Relaxation parameter.
   double relaxation_parameter = 1.0;
 
@@ -168,14 +135,11 @@ class NewtonSolver
   std::function<void(const NewtonSolver& solver, const Vec dx, Vec x)>
       _update_solution;
 
-  // Accumulated number of Krylov iterations since solve began
-  int _krylov_iterations;
-
   // Number of iterations
-  int _iteration;
+  // int _iteration;
 
   // Most recent residual and initial residual
-  double _residual, _residual0;
+  // double _residual, _residual0;
 
   // Linear solver
   dolfinx::la::petsc::KrylovSolver _solver;
diff --git a/cpp/demo/hyperelasticity/hyperelasticity.py b/cpp/demo/hyperelasticity/hyperelasticity.py
index 97d6d6ccd6e..995e755f97a 100644
--- a/cpp/demo/hyperelasticity/hyperelasticity.py
+++ b/cpp/demo/hyperelasticity/hyperelasticity.py
@@ -10,7 +10,6 @@
     Constant,
     FunctionSpace,
     Identity,
-    Measure,
     Mesh,
     TestFunction,
     TrialFunction,
@@ -18,6 +17,7 @@
     det,
     diff,
     ds,
+    dx,
     grad,
     inner,
     ln,
@@ -27,7 +27,7 @@
 
 # Function spaces
 e = element("Lagrange", "tetrahedron", 1, shape=(3,))
-e2 = element("Lagrange", "tetrahedron", 2, shape=(3,))
+e2 = element("Lagrange", "tetrahedron", 1, shape=(3,))
 mesh = Mesh(e)
 V = FunctionSpace(mesh, e2)
 
@@ -76,7 +76,7 @@
 psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
 
 
-dx = Measure("dx", metadata={"quadrature_degree": 3})
+# dx = Measure("dx", metadata={"quadrature_degree": 3})
 # Total potential energy
 Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds
 
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index a4d35edcc61..3c87c3133a6 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -162,7 +162,7 @@ int main(int argc, char* argv[])
         mesh::create_cell_partitioner(mesh::GhostMode::none)));
 
     auto element = basix::create_element<U>(
-        basix::element::family::P, basix::cell::type::tetrahedron, 2,
+        basix::element::family::P, basix::cell::type::tetrahedron, 1,
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 

From a1a1a516254f00af8aaf9b880b78d24a2c00284f Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 12:06:45 +0100
Subject: [PATCH 3/7] Incorporate solver into main.cpp

---
 cpp/demo/hyperelasticity/CMakeLists.txt |   2 +-
 cpp/demo/hyperelasticity/NewtonSolver.h |   6 -
 cpp/demo/hyperelasticity/main.cpp       | 231 +++++++++++++++++-------
 3 files changed, 171 insertions(+), 68 deletions(-)

diff --git a/cpp/demo/hyperelasticity/CMakeLists.txt b/cpp/demo/hyperelasticity/CMakeLists.txt
index 77273759b1f..007b04e24c5 100644
--- a/cpp/demo/hyperelasticity/CMakeLists.txt
+++ b/cpp/demo/hyperelasticity/CMakeLists.txt
@@ -35,7 +35,7 @@ else()
 
   set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-  add_executable(${PROJECT_NAME} main.cpp NewtonSolver.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
+  add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
   target_link_libraries(${PROJECT_NAME} dolfinx)
 
   # Set C++20 standard
diff --git a/cpp/demo/hyperelasticity/NewtonSolver.h b/cpp/demo/hyperelasticity/NewtonSolver.h
index 82d6ff2811b..ca5838bc64f 100644
--- a/cpp/demo/hyperelasticity/NewtonSolver.h
+++ b/cpp/demo/hyperelasticity/NewtonSolver.h
@@ -135,12 +135,6 @@ class NewtonSolver
   std::function<void(const NewtonSolver& solver, const Vec dx, Vec x)>
       _update_solution;
 
-  // Number of iterations
-  // int _iteration;
-
-  // Most recent residual and initial residual
-  // double _residual, _residual0;
-
   // Linear solver
   dolfinx::la::petsc::KrylovSolver _solver;
 
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 3c87c3133a6..d8c8c647ecc 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -44,12 +44,13 @@ class HyperElasticProblem
 {
 public:
   /// Constructor
-  HyperElasticProblem(fem::Form<T>& L, fem::Form<T>& J,
+  HyperElasticProblem(MPI_Comm comm, fem::Form<T>& L, fem::Form<T>& J,
                       const std::vector<fem::DirichletBC<T>>& bcs)
       : _l(L), _j(J), _bcs(bcs.begin(), bcs.end()),
-        _b(L.function_spaces()[0]->dofmap()->index_map,
-           L.function_spaces()[0]->dofmap()->index_map_bs()),
-        _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false))
+        _b_vec(L.function_spaces()[0]->dofmap()->index_map,
+               L.function_spaces()[0]->dofmap()->index_map_bs()),
+        _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false)),
+        _solver(comm), _dx(nullptr), _comm(comm)
   {
     auto map = L.function_spaces()[0]->dofmap()->index_map;
     const int bs = L.function_spaces()[0]->dofmap()->index_map_bs();
@@ -59,82 +60,194 @@ class HyperElasticProblem
     std::int64_t size_global = bs * map->size_global();
     VecCreateGhostBlockWithArray(map->comm(), bs, size_local, size_global,
                                  ghosts.size(), ghosts.data(),
-                                 _b.array().data(), &_b_petsc);
+                                 _b_vec.array().data(), &_b);
+
+    // Create linear solver. Default to LU.
+    _solver.set_options_prefix("nls_solve_");
+    la::petsc::options::set("nls_solve_ksp_type", "preonly");
+    la::petsc::options::set("nls_solve_pc_type", "lu");
+    _solver.set_from_options();
+
+    _matJ = _matA.mat();
   }
 
   /// Destructor
   virtual ~HyperElasticProblem()
   {
-    if (_b_petsc)
-      VecDestroy(&_b_petsc);
+    if (_b)
+      VecDestroy(&_b);
+    if (_dx)
+      VecDestroy(&_dx);
   }
 
-  /// @brief  Form
-  /// @return
-  auto form()
+  std::pair<int, bool> solve(Vec x)
   {
-    return [](Vec x)
+    // Reset iteration counts
+    int iteration = 0;
+    int krylov_iterations = 0;
+    double residual = -1;
+    double residual0 = 0;
+
+    auto converged
+        = [&iteration, &residual0, this](const Vec r) -> std::pair<double, bool>
     {
-      VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
-      VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
+      PetscReal residual = 0;
+      VecNorm(r, NORM_2, &residual);
+
+      // Relative residual
+      const double relative_residual = residual / residual0;
+
+      // Output iteration number and residual
+      if (dolfinx::MPI::rank(_comm.comm()) == 0)
+      {
+        spdlog::info("Newton iteration {}"
+                     ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
+                     iteration, residual, atol, relative_residual, rtol);
+      }
+
+      // Return true if convergence criterion is met
+      if (relative_residual < rtol or residual < atol)
+        return {residual, true};
+      else
+        return {residual, false};
     };
+
+    form(x);
+    assert(_b);
+    F(x, _b);
+
+    // Check convergence
+    bool newton_converged = false;
+    std::tie(residual, newton_converged) = converged(_b);
+
+    _solver.set_operators(_matJ, _matJ);
+
+    MatCreateVecs(_matJ, &_dx, nullptr);
+
+    // Start iterations
+    while (!newton_converged and iteration < max_it)
+    {
+      // Compute Jacobian
+      assert(_matJ);
+      J(x, _matJ);
+
+      // Perform linear solve and update total number of Krylov iterations
+      krylov_iterations += _solver.solve(_dx, _b);
+
+      // Update solution
+      VecAXPY(x, -relaxation_parameter, _dx);
+
+      // Increment iteration count
+      ++iteration;
+
+      // Compute F
+      form(x);
+      F(x, _b);
+
+      // Initialize _residual0
+      if (iteration == 1)
+      {
+        PetscReal _r = 0;
+        VecNorm(_dx, NORM_2, &_r);
+        residual0 = _r;
+      }
+
+      // Test for convergence
+      std::tie(residual, newton_converged) = converged(_b);
+    }
+
+    if (newton_converged)
+    {
+      if (dolfinx::MPI::rank(_comm.comm()) == 0)
+      {
+        spdlog::info("Newton solver finished in {} iterations and {} linear "
+                     "solver iterations.",
+                     iteration, krylov_iterations);
+      }
+    }
+    else
+    {
+      throw std::runtime_error("Newton solver did not converge.");
+    }
+
+    return {iteration, newton_converged};
+  }
+
+  /// @brief  Form
+  /// @return
+  void form(Vec x)
+  {
+    VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
+    VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
   }
 
   /// Compute F at current point x
-  auto F()
+  void F(const Vec x, Vec)
   {
-    return [&](const Vec x, Vec)
-    {
-      // Assemble b and update ghosts
-      std::span b(_b.array());
-      std::ranges::fill(b, 0);
-      fem::assemble_vector(b, _l);
-      VecGhostUpdateBegin(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
-      VecGhostUpdateEnd(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
-
-      // Set bcs
-      Vec x_local;
-      VecGhostGetLocalForm(x, &x_local);
-      PetscInt n = 0;
-      VecGetSize(x_local, &n);
-      const T* _x = nullptr;
-      VecGetArrayRead(x_local, &_x);
-      std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
-                            { bc.get().set(b, x, -1); });
-      VecRestoreArrayRead(x_local, &_x);
-    };
+    // Assemble b and update ghosts
+    std::span b(_b_vec.array());
+    std::ranges::fill(b, 0);
+    fem::assemble_vector(b, _l);
+    VecGhostUpdateBegin(_b, ADD_VALUES, SCATTER_REVERSE);
+    VecGhostUpdateEnd(_b, ADD_VALUES, SCATTER_REVERSE);
+
+    // Set bcs
+    Vec x_local;
+    VecGhostGetLocalForm(x, &x_local);
+    PetscInt n = 0;
+    VecGetSize(x_local, &n);
+    const T* _x = nullptr;
+    VecGetArrayRead(x_local, &_x);
+    std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
+                          { bc.get().set(b, x, -1); });
+    VecRestoreArrayRead(x_local, &_x);
   }
 
   /// Compute J = F' at current point x
-  auto J()
+  void J(const Vec, Mat A)
   {
-    return [&](const Vec, Mat A)
-    {
-      MatZeroEntries(A);
-      fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), _j,
-                           _bcs);
-      MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
-      MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
-      fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
-                        *_j.function_spaces()[0], _bcs);
-      MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
-      MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
-    };
+    MatZeroEntries(A);
+    fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), _j,
+                         _bcs);
+    MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
+    MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
+    fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
+                      *_j.function_spaces()[0], _bcs);
+    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
   }
 
-  /// RHS vector
-  Vec vector() { return _b_petsc; }
+  /// @brief Relative convergence tolerance.
+  double rtol = 1e-9;
 
-  /// Jacobian matrix
-  Mat matrix() { return _matA.mat(); }
+  /// @brief Absolute convergence tolerance.
+  double atol = 1e-10;
 
 private:
   fem::Form<T>& _l;
   fem::Form<T>& _j;
   std::vector<std::reference_wrapper<const fem::DirichletBC<T>>> _bcs;
-  la::Vector<T> _b;
-  Vec _b_petsc = nullptr;
+  la::Vector<T> _b_vec;
+  Vec _b = nullptr;
+
+  // Jacobian matrix
   la::petsc::Matrix _matA;
+  Mat _matJ = nullptr;
+
+  // Linear solver
+  dolfinx::la::petsc::KrylovSolver _solver;
+
+  // Solution vector
+  Vec _dx = nullptr;
+
+  // MPI communicator
+  dolfinx::MPI::Comm _comm;
+
+  /// @brief Maximum number of iterations.
+  int max_it = 50;
+
+  /// @brief Relaxation parameter.
+  double relaxation_parameter = 1.0;
 };
 
 int main(int argc, char* argv[])
@@ -246,16 +359,12 @@ int main(int argc, char* argv[])
         = {fem::DirichletBC<T>(std::vector<T>{0, 0, 0}, bdofs_left, V),
            fem::DirichletBC<T>(u_rotation, bdofs_right)};
 
-    HyperElasticProblem problem(L, a, bcs);
-    NewtonSolver newton_solver(mesh->comm());
-    newton_solver.setF(problem.F(), problem.vector());
-    newton_solver.setJ(problem.J(), problem.matrix());
-    newton_solver.set_form(problem.form());
-    newton_solver.rtol = 10 * std::numeric_limits<T>::epsilon();
-    newton_solver.atol = 10 * std::numeric_limits<T>::epsilon();
+    HyperElasticProblem problem(mesh->comm(), L, a, bcs);
+    problem.rtol = 10 * std::numeric_limits<T>::epsilon();
+    problem.atol = 10 * std::numeric_limits<T>::epsilon();
 
     la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
-    auto [niter, success] = newton_solver.solve(_u.vec());
+    auto [niter, success] = problem.solve(_u.vec());
     std::cout << "Number of Newton iterations: " << niter << std::endl;
 
     // Compute Cauchy stress. Construct appropriate Basix element for

From bf2fefc2ba1af5354d71e5287e892f348a987e4a Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 12:36:48 +0100
Subject: [PATCH 4/7] Reduce class variables

---
 cpp/demo/hyperelasticity/main.cpp | 26 ++++++++++----------------
 1 file changed, 10 insertions(+), 16 deletions(-)

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index d8c8c647ecc..fb888c9ddb7 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -50,7 +50,7 @@ class HyperElasticProblem
         _b_vec(L.function_spaces()[0]->dofmap()->index_map,
                L.function_spaces()[0]->dofmap()->index_map_bs()),
         _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false)),
-        _solver(comm), _dx(nullptr), _comm(comm)
+        _solver(comm), _comm(comm)
   {
     auto map = L.function_spaces()[0]->dofmap()->index_map;
     const int bs = L.function_spaces()[0]->dofmap()->index_map_bs();
@@ -76,8 +76,6 @@ class HyperElasticProblem
   {
     if (_b)
       VecDestroy(&_b);
-    if (_dx)
-      VecDestroy(&_dx);
   }
 
   std::pair<int, bool> solve(Vec x)
@@ -122,8 +120,10 @@ class HyperElasticProblem
 
     _solver.set_operators(_matJ, _matJ);
 
-    MatCreateVecs(_matJ, &_dx, nullptr);
+    Vec dx;
+    MatCreateVecs(_matJ, &dx, nullptr);
 
+    int max_it = 50;
     // Start iterations
     while (!newton_converged and iteration < max_it)
     {
@@ -132,10 +132,11 @@ class HyperElasticProblem
       J(x, _matJ);
 
       // Perform linear solve and update total number of Krylov iterations
-      krylov_iterations += _solver.solve(_dx, _b);
+      krylov_iterations += _solver.solve(dx, _b);
 
       // Update solution
-      VecAXPY(x, -relaxation_parameter, _dx);
+      double relaxation_parameter = 1.0;
+      VecAXPY(x, -relaxation_parameter, dx);
 
       // Increment iteration count
       ++iteration;
@@ -148,7 +149,7 @@ class HyperElasticProblem
       if (iteration == 1)
       {
         PetscReal _r = 0;
-        VecNorm(_dx, NORM_2, &_r);
+        VecNorm(dx, NORM_2, &_r);
         residual0 = _r;
       }
 
@@ -170,6 +171,8 @@ class HyperElasticProblem
       throw std::runtime_error("Newton solver did not converge.");
     }
 
+    VecDestroy(&dx);
+
     return {iteration, newton_converged};
   }
 
@@ -237,17 +240,8 @@ class HyperElasticProblem
   // Linear solver
   dolfinx::la::petsc::KrylovSolver _solver;
 
-  // Solution vector
-  Vec _dx = nullptr;
-
   // MPI communicator
   dolfinx::MPI::Comm _comm;
-
-  /// @brief Maximum number of iterations.
-  int max_it = 50;
-
-  /// @brief Relaxation parameter.
-  double relaxation_parameter = 1.0;
 };
 
 int main(int argc, char* argv[])

From 10388364d203db9f5607e8930fffe1d4c0df5d30 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 13:53:03 +0100
Subject: [PATCH 5/7] Simplify

---
 cpp/demo/hyperelasticity/main.cpp | 48 ++++++++++++-------------------
 1 file changed, 18 insertions(+), 30 deletions(-)

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index fb888c9ddb7..32827271303 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -44,13 +44,13 @@ class HyperElasticProblem
 {
 public:
   /// Constructor
-  HyperElasticProblem(MPI_Comm comm, fem::Form<T>& L, fem::Form<T>& J,
+  HyperElasticProblem(fem::Form<T>& L, fem::Form<T>& J,
                       const std::vector<fem::DirichletBC<T>>& bcs)
       : _l(L), _j(J), _bcs(bcs.begin(), bcs.end()),
         _b_vec(L.function_spaces()[0]->dofmap()->index_map,
                L.function_spaces()[0]->dofmap()->index_map_bs()),
-        _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false)),
-        _solver(comm), _comm(comm)
+        _matJ(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false)),
+        _solver(L.function_spaces()[0]->dofmap()->index_map->comm())
   {
     auto map = L.function_spaces()[0]->dofmap()->index_map;
     const int bs = L.function_spaces()[0]->dofmap()->index_map_bs();
@@ -67,8 +67,6 @@ class HyperElasticProblem
     la::petsc::options::set("nls_solve_ksp_type", "preonly");
     la::petsc::options::set("nls_solve_pc_type", "lu");
     _solver.set_from_options();
-
-    _matJ = _matA.mat();
   }
 
   /// Destructor
@@ -96,12 +94,9 @@ class HyperElasticProblem
       const double relative_residual = residual / residual0;
 
       // Output iteration number and residual
-      if (dolfinx::MPI::rank(_comm.comm()) == 0)
-      {
-        spdlog::info("Newton iteration {}"
-                     ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
-                     iteration, residual, atol, relative_residual, rtol);
-      }
+      spdlog::info("Newton iteration {}"
+                   ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
+                   iteration, residual, atol, relative_residual, rtol);
 
       // Return true if convergence criterion is met
       if (relative_residual < rtol or residual < atol)
@@ -112,24 +107,24 @@ class HyperElasticProblem
 
     form(x);
     assert(_b);
-    F(x, _b);
+    F(x);
 
     // Check convergence
     bool newton_converged = false;
     std::tie(residual, newton_converged) = converged(_b);
 
-    _solver.set_operators(_matJ, _matJ);
+    _solver.set_operators(_matJ.mat(), _matJ.mat());
 
     Vec dx;
-    MatCreateVecs(_matJ, &dx, nullptr);
+    MatCreateVecs(_matJ.mat(), &dx, nullptr);
 
     int max_it = 50;
     // Start iterations
     while (!newton_converged and iteration < max_it)
     {
       // Compute Jacobian
-      assert(_matJ);
-      J(x, _matJ);
+      assert(_matJ.mat());
+      J(x, _matJ.mat());
 
       // Perform linear solve and update total number of Krylov iterations
       krylov_iterations += _solver.solve(dx, _b);
@@ -143,7 +138,7 @@ class HyperElasticProblem
 
       // Compute F
       form(x);
-      F(x, _b);
+      F(x);
 
       // Initialize _residual0
       if (iteration == 1)
@@ -159,12 +154,9 @@ class HyperElasticProblem
 
     if (newton_converged)
     {
-      if (dolfinx::MPI::rank(_comm.comm()) == 0)
-      {
-        spdlog::info("Newton solver finished in {} iterations and {} linear "
-                     "solver iterations.",
-                     iteration, krylov_iterations);
-      }
+      spdlog::info("Newton solver finished in {} iterations and {} linear "
+                   "solver iterations.",
+                   iteration, krylov_iterations);
     }
     else
     {
@@ -185,7 +177,7 @@ class HyperElasticProblem
   }
 
   /// Compute F at current point x
-  void F(const Vec x, Vec)
+  void F(const Vec x)
   {
     // Assemble b and update ghosts
     std::span b(_b_vec.array());
@@ -234,14 +226,10 @@ class HyperElasticProblem
   Vec _b = nullptr;
 
   // Jacobian matrix
-  la::petsc::Matrix _matA;
-  Mat _matJ = nullptr;
+  la::petsc::Matrix _matJ;
 
   // Linear solver
   dolfinx::la::petsc::KrylovSolver _solver;
-
-  // MPI communicator
-  dolfinx::MPI::Comm _comm;
 };
 
 int main(int argc, char* argv[])
@@ -353,7 +341,7 @@ int main(int argc, char* argv[])
         = {fem::DirichletBC<T>(std::vector<T>{0, 0, 0}, bdofs_left, V),
            fem::DirichletBC<T>(u_rotation, bdofs_right)};
 
-    HyperElasticProblem problem(mesh->comm(), L, a, bcs);
+    HyperElasticProblem problem(L, a, bcs);
     problem.rtol = 10 * std::numeric_limits<T>::epsilon();
     problem.atol = 10 * std::numeric_limits<T>::epsilon();
 

From 877f5831d4a855eec18da2d183cbc8f8270d42b4 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 14:15:27 +0100
Subject: [PATCH 6/7] Simplify

---
 cpp/demo/hyperelasticity/main.cpp | 25 +++++++++----------------
 1 file changed, 9 insertions(+), 16 deletions(-)

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 32827271303..1fd441ea368 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -76,12 +76,12 @@ class HyperElasticProblem
       VecDestroy(&_b);
   }
 
+  /// @brief Newton Solver
+  /// @param x Solution vector
+  /// @return Iteration count and flag indicating convergence
   std::pair<int, bool> solve(Vec x)
   {
-    // Reset iteration counts
     int iteration = 0;
-    int krylov_iterations = 0;
-    double residual = -1;
     double residual0 = 0;
 
     auto converged
@@ -105,13 +105,10 @@ class HyperElasticProblem
         return {residual, false};
     };
 
-    form(x);
     assert(_b);
     F(x);
 
-    // Check convergence
-    bool newton_converged = false;
-    std::tie(residual, newton_converged) = converged(_b);
+    auto [residual, newton_converged] = converged(_b);
 
     _solver.set_operators(_matJ.mat(), _matJ.mat());
 
@@ -119,6 +116,8 @@ class HyperElasticProblem
     MatCreateVecs(_matJ.mat(), &dx, nullptr);
 
     int max_it = 50;
+    int krylov_iterations = 0;
+
     // Start iterations
     while (!newton_converged and iteration < max_it)
     {
@@ -137,10 +136,9 @@ class HyperElasticProblem
       ++iteration;
 
       // Compute F
-      form(x);
       F(x);
 
-      // Initialize _residual0
+      // Initialize residual0
       if (iteration == 1)
       {
         PetscReal _r = 0;
@@ -168,17 +166,12 @@ class HyperElasticProblem
     return {iteration, newton_converged};
   }
 
-  /// @brief  Form
-  /// @return
-  void form(Vec x)
+  /// Compute F at current point x
+  void F(const Vec x)
   {
     VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
     VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
-  }
 
-  /// Compute F at current point x
-  void F(const Vec x)
-  {
     // Assemble b and update ghosts
     std::span b(_b_vec.array());
     std::ranges::fill(b, 0);

From 83cb6c3bc2dc55d69d9b05f66c6e425fd6467305 Mon Sep 17 00:00:00 2001
From: Chris Richardson <chris@bpi.cam.ac.uk>
Date: Wed, 17 Sep 2025 14:49:32 +0100
Subject: [PATCH 7/7] Remove external files

---
 cpp/demo/hyperelasticity/NewtonSolver.cpp   | 168 --------------------
 cpp/demo/hyperelasticity/NewtonSolver.h     | 148 -----------------
 cpp/demo/hyperelasticity/hyperelasticity.py |   5 +-
 cpp/demo/hyperelasticity/main.cpp           |  15 +-
 4 files changed, 5 insertions(+), 331 deletions(-)
 delete mode 100644 cpp/demo/hyperelasticity/NewtonSolver.cpp
 delete mode 100644 cpp/demo/hyperelasticity/NewtonSolver.h

diff --git a/cpp/demo/hyperelasticity/NewtonSolver.cpp b/cpp/demo/hyperelasticity/NewtonSolver.cpp
deleted file mode 100644
index 3e76388d285..00000000000
--- a/cpp/demo/hyperelasticity/NewtonSolver.cpp
+++ /dev/null
@@ -1,168 +0,0 @@
-// Copyright (C) 2005-2021 Garth N. Wells
-//
-// This file is part of DOLFINx (https://www.fenicsproject.org)
-//
-// SPDX-License-Identifier:    LGPL-3.0-or-later
-
-// #ifdef HAS_PETSC
-
-#include "NewtonSolver.h"
-#include <dolfinx/common/MPI.h>
-#include <dolfinx/common/log.h>
-#include <dolfinx/la/petsc.h>
-#include <string>
-#include <utility>
-
-using namespace dolfinx;
-
-//-----------------------------------------------------------------------------
-NewtonSolver::NewtonSolver(MPI_Comm comm)
-    : _solver(comm), _dx(nullptr), _comm(comm)
-{
-  // Create linear solver if not already created. Default to LU.
-  _solver.set_options_prefix("nls_solve_");
-  la::petsc::options::set("nls_solve_ksp_type", "preonly");
-  la::petsc::options::set("nls_solve_pc_type", "lu");
-  _solver.set_from_options();
-}
-//-----------------------------------------------------------------------------
-NewtonSolver::~NewtonSolver()
-{
-  if (_b)
-    VecDestroy(&_b);
-  if (_dx)
-    VecDestroy(&_dx);
-  if (_matJ)
-    MatDestroy(&_matJ);
-  if (_matP)
-    MatDestroy(&_matP);
-}
-//-----------------------------------------------------------------------------
-void NewtonSolver::setF(std::function<void(const Vec, Vec)> F,
-                                    Vec b)
-{
-  _fnF = std::move(F);
-  _b = b;
-  PetscObjectReference((PetscObject)_b);
-}
-//-----------------------------------------------------------------------------
-void NewtonSolver::setJ(std::function<void(const Vec, Mat)> J,
-                                    Mat Jmat)
-{
-  _fnJ = std::move(J);
-  _matJ = Jmat;
-  PetscObjectReference((PetscObject)_matJ);
-}
-//-----------------------------------------------------------------------------
-void NewtonSolver::set_form(std::function<void(Vec)> form)
-{
-  _system = std::move(form);
-}
-//-----------------------------------------------------------------------------
-std::pair<int, bool> NewtonSolver::solve(Vec x)
-{
-  // Reset iteration counts
-  int iteration = 0;
-  int krylov_iterations = 0;
-  double residual = -1;
-  double residual0 = 0;
-
-  auto converged = [&iteration, &residual, &residual0,
-                    this](const Vec r) -> std::pair<double, bool>
-  {
-    PetscReal residual = 0;
-    VecNorm(r, NORM_2, &residual);
-
-    // Relative residual
-    const double relative_residual = residual / residual0;
-
-    // Output iteration number and residual
-    if (dolfinx::MPI::rank(_comm.comm()) == 0)
-    {
-      spdlog::info("Newton iteration {}"
-                   ": r (abs) = {} (tol = {}), r (rel) = {} (tol = {})",
-                   iteration, residual, atol, relative_residual, rtol);
-    }
-
-    // Return true if convergence criterion is met
-    if (relative_residual < rtol or residual < atol)
-      return {residual, true};
-    else
-      return {residual, false};
-  };
-
-  if (!_fnF)
-  {
-    throw std::runtime_error("Function for computing residual vector has not "
-                             "been provided to the NewtonSolver.");
-  }
-
-  if (!_fnJ)
-  {
-    throw std::runtime_error("Function for computing Jacobian has not "
-                             "been provided to the NewtonSolver.");
-  }
-
-  if (_system)
-    _system(x);
-  assert(_b);
-  _fnF(x, _b);
-
-  // Check convergence
-  bool newton_converged = false;
-  std::tie(residual, newton_converged) = converged(_b);
-
-  _solver.set_operators(_matJ, _matJ);
-
-  MatCreateVecs(_matJ, &_dx, nullptr);
-
-  // Start iterations
-  while (!newton_converged and iteration < max_it)
-  {
-    // Compute Jacobian
-    assert(_matJ);
-    _fnJ(x, _matJ);
-
-    // Perform linear solve and update total number of Krylov iterations
-    krylov_iterations += _solver.solve(_dx, _b);
-
-    // Update solution
-    VecAXPY(x, -relaxation_parameter, _dx);
-
-    // Increment iteration count
-    ++iteration;
-
-    // Compute F
-    if (_system)
-      _system(x);
-    _fnF(x, _b);
-
-    // Initialize _residual0
-    if (iteration == 1)
-    {
-      PetscReal _r = 0;
-      VecNorm(_dx, NORM_2, &_r);
-      residual0 = _r;
-    }
-
-    // Test for convergence
-    std::tie(residual, newton_converged) = converged(_b);
-  }
-
-  if (newton_converged)
-  {
-    if (dolfinx::MPI::rank(_comm.comm()) == 0)
-    {
-      spdlog::info("Newton solver finished in {} iterations and {} linear "
-                   "solver iterations.",
-                   iteration, krylov_iterations);
-    }
-  }
-  else
-  {
-    throw std::runtime_error("Newton solver did not converge.");
-  }
-
-  return {iteration, newton_converged};
-}
-// #endif
diff --git a/cpp/demo/hyperelasticity/NewtonSolver.h b/cpp/demo/hyperelasticity/NewtonSolver.h
deleted file mode 100644
index ca5838bc64f..00000000000
--- a/cpp/demo/hyperelasticity/NewtonSolver.h
+++ /dev/null
@@ -1,148 +0,0 @@
-// Copyright (C) 2005-2021 Garth N. Wells
-//
-// This file is part of DOLFINx (https://www.fenicsproject.org)
-//
-// SPDX-License-Identifier:    LGPL-3.0-or-later
-
-#pragma once
-
-// #ifdef HAS_PETSC
-
-#include <dolfinx/common/MPI.h>
-#include <dolfinx/la/petsc.h>
-#include <functional>
-#include <memory>
-#include <petscmat.h>
-#include <petscvec.h>
-#include <utility>
-
-
-/// @brief A Newton solver for nonlinear systems of equations of the
-/// form \f$F(x) = 0\f$.
-///
-/// It solves \f[ \left. \frac{dF}{dx} \right|_{x} \Delta x = F(x) \f]
-/// with default update \f$x \leftarrow x - \Delta x\f$.
-///
-/// It relies on PETSc for linear algebra backends.
-class NewtonSolver
-{
-public:
-  /// @brief Create nonlinear solver
-  /// @param[in] comm MPI communicator for the solver
-  explicit NewtonSolver(MPI_Comm comm);
-
-  /// Move constructor
-  NewtonSolver(NewtonSolver&& solver) = default;
-
-  // Copy constructor (deleted)
-  NewtonSolver(const NewtonSolver& solver) = delete;
-
-  /// Move assignment constructor
-  NewtonSolver& operator=(NewtonSolver&& solver) = default;
-
-  // Assignment operator (deleted)
-  NewtonSolver& operator=(const NewtonSolver& solver) = delete;
-
-  /// Destructor
-  ~NewtonSolver();
-
-  /// @brief Set the function for computing the residual \f$F(x) = 0\f$
-  /// and the vector to assemble the residual into.
-  /// @param[in] F Function to compute/assemble the residual vector `b`.
-  /// The first argument to the function is the solution vector `x` and
-  /// the second is the vector `b` to assemble into.
-  /// @param[in] b Vector to assemble to residual into.
-  void setF(std::function<void(const Vec, Vec)> F, Vec b);
-
-  /// @brief Set the function for computing the Jacobian \f$J:=dF/dx\f$
-  /// and the matrix to assemble the Jacobian into.
-  /// @param[in] J Function to compute the Jacobian matrix `Jmat`. The
-  /// first argument to the function is the solution vector `x` and the
-  /// second is the matrix to assemble into.
-  /// @param[in] Jmat Matrix to assemble the Jacobian into.
-  void setJ(std::function<void(const Vec, Mat)> J, Mat Jmat);
-
-  /// @brief Set the function that is called before the residual or
-  /// Jacobian are computed. It is commonly used to update ghost values.
-  ///
-  /// @param[in] form Function to call. It takes the (latest) solution
-  /// vector `x` as an argument.
-  void set_form(std::function<void(Vec)> form);
-
-  /// @brief Optional set function that is called after each inner solve for
-  /// the Newton increment to update the solution.
-  ///
-  /// The function `update` takes `this`, the Newton increment `dx`, and
-  /// the vector `x` from the start of the Newton iteration.
-  ///
-  /// By default, the update is x <- x - dx
-  ///
-  /// @param[in] update The function that updates the solution
-  void set_update(
-      std::function<void(const NewtonSolver& solver, const Vec, Vec)> update);
-
-  /// @brief Solve the nonlinear problem.
-  ///
-  /// @param[in,out] x The solution vector. It should be set the initial
-  /// solution guess.
-  /// @return (number of Newton iterations, whether iteration converged)
-  std::pair<int, bool> solve(Vec x);
-
-  /// @brief Maximum number of iterations.
-  int max_it = 50;
-
-  /// @brief Relative convergence tolerance.
-  double rtol = 1e-9;
-
-  /// @brief Absolute convergence tolerance.
-  double atol = 1e-10;
-
-  /// @brief Relaxation parameter.
-  double relaxation_parameter = 1.0;
-
-private:
-  // Function for computing the residual vector. The first argument is
-  // the latest solution vector x and the second argument is the
-  // residual vector.
-  std::function<void(const Vec x, Vec b)> _fnF;
-
-  // Function for computing the Jacobian matrix operator. The first
-  // argument is the latest solution vector x and the second argument is
-  // the matrix operator.
-  std::function<void(const Vec x, Mat J)> _fnJ;
-
-  // Function for computing the preconditioner matrix operator. The
-  // first argument is the latest solution vector x and the second
-  // argument is the matrix operator.
-  std::function<void(const Vec x, Mat P)> _fnP;
-
-  // Function called before the residual and Jacobian function at each
-  // iteration.
-  std::function<void(const Vec x)> _system;
-
-  // Residual vector
-  Vec _b = nullptr;
-
-  // Jacobian matrix and preconditioner matrix
-  Mat _matJ = nullptr, _matP = nullptr;
-
-  // Function to check for convergence
-  std::function<std::pair<double, bool>(const NewtonSolver& solver,
-                                        const Vec r)>
-      _converged;
-
-  // Function to update the solution after inner solve for the Newton increment
-  std::function<void(const NewtonSolver& solver, const Vec dx, Vec x)>
-      _update_solution;
-
-  // Linear solver
-  dolfinx::la::petsc::KrylovSolver _solver;
-
-  // Solution vector
-  Vec _dx = nullptr;
-
-  // MPI communicator
-  dolfinx::MPI::Comm _comm;
-};
-
-// #endif
diff --git a/cpp/demo/hyperelasticity/hyperelasticity.py b/cpp/demo/hyperelasticity/hyperelasticity.py
index 995e755f97a..f361c4171ab 100644
--- a/cpp/demo/hyperelasticity/hyperelasticity.py
+++ b/cpp/demo/hyperelasticity/hyperelasticity.py
@@ -27,9 +27,8 @@
 
 # Function spaces
 e = element("Lagrange", "tetrahedron", 1, shape=(3,))
-e2 = element("Lagrange", "tetrahedron", 1, shape=(3,))
 mesh = Mesh(e)
-V = FunctionSpace(mesh, e2)
+V = FunctionSpace(mesh, e)
 
 # Trial and test functions
 du = TrialFunction(V)  # Incremental displacement
@@ -75,8 +74,6 @@
 # Stored strain energy density (compressible neo-Hookean model)
 psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
 
-
-# dx = Measure("dx", metadata={"quadrature_degree": 3})
 # Total potential energy
 Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds
 
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 1fd441ea368..83ea9be399c 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -14,7 +14,6 @@
 
 // ## C++ program
 
-#include "NewtonSolver.h"
 #include "hyperelasticity.h"
 #include <algorithm>
 #include <basix/finite-element.h>
@@ -72,8 +71,8 @@ class HyperElasticProblem
   /// Destructor
   virtual ~HyperElasticProblem()
   {
-    if (_b)
-      VecDestroy(&_b);
+    assert(_b);
+    VecDestroy(&_b);
   }
 
   /// @brief Newton Solver
@@ -82,7 +81,7 @@ class HyperElasticProblem
   std::pair<int, bool> solve(Vec x)
   {
     int iteration = 0;
-    double residual0 = 0;
+    PetscReal residual0 = 0;
 
     auto converged
         = [&iteration, &residual0, this](const Vec r) -> std::pair<double, bool>
@@ -140,11 +139,7 @@ class HyperElasticProblem
 
       // Initialize residual0
       if (iteration == 1)
-      {
-        PetscReal _r = 0;
-        VecNorm(dx, NORM_2, &_r);
-        residual0 = _r;
-      }
+        VecNorm(dx, NORM_2, &residual0);
 
       // Test for convergence
       std::tie(residual, newton_converged) = converged(_b);
@@ -157,9 +152,7 @@ class HyperElasticProblem
                    iteration, krylov_iterations);
     }
     else
-    {
       throw std::runtime_error("Newton solver did not converge.");
-    }
 
     VecDestroy(&dx);