From e41b0b15cf059122119bf850460a8b08527e5fa7 Mon Sep 17 00:00:00 2001
From: "Garth N. Wells" <gnw20@cam.ac.uk>
Date: Thu, 26 Sep 2024 15:49:32 +0100
Subject: [PATCH] Use `std::optional` to reduce code duplication in boundary
 condition setting (#3434)

* Use std::optional in bc setting

* Small update

* Demo update

* Binding fix

* Add no convert

* Ruff fix

* Reduce code duplication

* Small fix

* Tidy up

* Tidy

* Named argument update

* Test updates

* Small fix

* More fixes

* Add BC set function

* Simplify

* Updates

* C++ updates

* Remove commented code

* Updates

* Fix warning

* Demo update

* Improve docs

* Small updates

* Simplifications

* Argument name fix in deprecated function

* Add Python docs

* Remove duplicated code

* Simplification
---
 cpp/demo/biharmonic/main.cpp                  |   2 +-
 cpp/demo/codim_0_assembly/main.cpp            |   6 +-
 cpp/demo/hyperelasticity/main.cpp             |  13 +-
 .../interpolation_different_meshes/main.cpp   |   3 +-
 cpp/demo/poisson/main.cpp                     |   2 +-
 cpp/demo/poisson_matrix_free/main.cpp         |   8 +-
 cpp/dolfinx/fem/DirichletBC.h                 | 202 +++++++++---------
 cpp/dolfinx/fem/assemble_vector_impl.h        |  85 ++++----
 cpp/dolfinx/fem/assembler.h                   |  61 +-----
 cpp/dolfinx/fem/petsc.h                       |  34 +--
 python/demo/demo_elasticity.py                |   4 +-
 python/demo/demo_poisson_matrix_free.py       |   8 +-
 python/demo/demo_pyamg.py                     |   5 +-
 python/demo/demo_static-condensation.py       |   4 +-
 python/demo/demo_stokes.py                    |   3 +-
 python/demo/demo_types.py                     |   4 +-
 python/dolfinx/fem/__init__.py                |   2 +-
 python/dolfinx/fem/assemble.py                |  21 +-
 python/dolfinx/fem/bcs.py                     |  35 ++-
 python/dolfinx/fem/petsc.py                   |  56 ++---
 python/dolfinx/wrappers/assemble.cpp          |  30 +--
 python/dolfinx/wrappers/fem.cpp               |  22 +-
 python/test/unit/fem/test_assemble_domains.py |   8 +-
 python/test/unit/fem/test_assemble_submesh.py |   2 +-
 python/test/unit/fem/test_bcs.py              |  25 +--
 .../unit/fem/test_custom_basix_element.py     |   3 +-
 python/test/unit/fem/test_fem_pipeline.py     |   6 +-
 .../fem/test_petsc_nonlinear_assembler.py     |  16 +-
 python/test/unit/nls/test_newton.py           |   4 +-
 29 files changed, 334 insertions(+), 340 deletions(-)

diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
index a9a61d8631f..4dec92293bc 100644
--- a/cpp/demo/biharmonic/main.cpp
+++ b/cpp/demo/biharmonic/main.cpp
@@ -239,7 +239,7 @@ int main(int argc, char* argv[])
     fem::assemble_vector(b.mutable_array(), *L);
     fem::apply_lifting<T, U>(b.mutable_array(), {a}, {{bc}}, {}, T(1.0));
     b.scatter_rev(std::plus<T>());
-    fem::set_bc<T, U>(b.mutable_array(), {bc});
+    bc->set(b.mutable_array(), std::nullopt);
 
     la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
     la::petsc::options::set("ksp_type", "preonly");
diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp
index 29c91a21a92..95537e79b29 100644
--- a/cpp/demo/codim_0_assembly/main.cpp
+++ b/cpp/demo/codim_0_assembly/main.cpp
@@ -101,7 +101,7 @@ int main(int argc, char* argv[])
     std::map<std::shared_ptr<const mesh::Mesh<U>>,
              std::span<const std::int32_t>>
         entity_maps
-        = {{const_ptr, std::span<const std::int32_t>(mesh_to_submesh.data(),
+        = {{const_ptr, std::span(mesh_to_submesh.data(),
                                                      mesh_to_submesh.size())}};
 
     // Next we compute the integration entities on the integration domain `mesh`
@@ -113,8 +113,8 @@ int main(int argc, char* argv[])
         fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2), tdim);
 
     subdomain_map[fem::IntegralType::cell].push_back(
-        {3, std::span<const std::int32_t>(integration_entities.data(),
-                                          integration_entities.size())});
+        {3,
+         std::span(integration_entities.data(), integration_entities.size())});
 
     // We can now create the bi-linear form
     auto a_mixed = std::make_shared<fem::Form<T>>(
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
index 9734dd640d6..01790d7408e 100644
--- a/cpp/demo/hyperelasticity/main.cpp
+++ b/cpp/demo/hyperelasticity/main.cpp
@@ -86,8 +86,8 @@ class HyperElasticProblem
     return [&](const Vec x, Vec)
     {
       // Assemble b and update ghosts
-      std::span<T> b(_b.mutable_array());
-      std::ranges::fill(b, 0.0);
+      std::span b(_b.mutable_array());
+      std::ranges::fill(b, 0);
       fem::assemble_vector<T>(b, *_l);
       VecGhostUpdateBegin(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
       VecGhostUpdateEnd(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
@@ -97,10 +97,11 @@ class HyperElasticProblem
       VecGhostGetLocalForm(x, &x_local);
       PetscInt n = 0;
       VecGetSize(x_local, &n);
-      const T* array = nullptr;
-      VecGetArrayRead(x_local, &array);
-      fem::set_bc<T>(b, _bcs, std::span<const T>(array, n), -1.0);
-      VecRestoreArrayRead(x, &array);
+      const T* _x = nullptr;
+      VecGetArrayRead(x_local, &_x);
+      std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
+                            { bc->set(b, x, -1); });
+      VecRestoreArrayRead(x_local, &_x);
     };
   }
 
diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
index c7d381f4c3c..c83cd0852fc 100644
--- a/cpp/demo/interpolation_different_meshes/main.cpp
+++ b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -78,8 +78,7 @@ int main(int argc, char* argv[])
         = fem::create_interpolation_data(
             u_hex->function_space()->mesh()->geometry(),
             *u_hex->function_space()->element(),
-            *u_tet->function_space()->mesh(),
-            std::span<const std::int32_t>(cells), 1e-8);
+            *u_tet->function_space()->mesh(), std::span(cells), 1e-8);
     u_hex->interpolate(*u_tet, cells, interpolation_data);
 
 #ifdef HAS_ADIOS2
diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
index 01f14418285..6d5d6624d4e 100644
--- a/cpp/demo/poisson/main.cpp
+++ b/cpp/demo/poisson/main.cpp
@@ -225,7 +225,7 @@ int main(int argc, char* argv[])
     fem::assemble_vector(b.mutable_array(), *L);
     fem::apply_lifting<T, U>(b.mutable_array(), {a}, {{bc}}, {}, T(1));
     b.scatter_rev(std::plus<T>());
-    fem::set_bc<T, U>(b.mutable_array(), {bc});
+    bc->set(b.mutable_array(), std::nullopt);
 
     la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
     la::petsc::options::set("ksp_type", "preonly");
diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp
index 31d2c6fc531..df7bb19aa41 100644
--- a/cpp/demo/poisson_matrix_free/main.cpp
+++ b/cpp/demo/poisson_matrix_free/main.cpp
@@ -182,14 +182,14 @@ void solver(MPI_Comm comm)
 
   // Apply lifting to account for Dirichlet boundary condition
   // b <- b - A * x_bc
-  fem::set_bc<T, U>(ui->x()->mutable_array(), {bc}, T(-1));
+  bc->set(ui->x()->mutable_array(), std::nullopt, T(-1));
   fem::assemble_vector(b.mutable_array(), *M);
 
   // Communicate ghost values
   b.scatter_rev(std::plus<T>());
 
   // Set BC dofs to zero (effectively zeroes columns of A)
-  fem::set_bc<T, U>(b.mutable_array(), {bc}, T(0));
+  bc->set(b.mutable_array(), std::nullopt, T(0));
 
   b.scatter_fwd();
 
@@ -212,7 +212,7 @@ void solver(MPI_Comm comm)
                          fem::make_coefficients_span(coeff));
 
     // Set BC dofs to zero (effectively zeroes rows of A)
-    fem::set_bc<T, U>(y.mutable_array(), {bc}, T(0));
+    bc->set(y.mutable_array(), std::nullopt, T(0));
 
     // Accumulate ghost values
     y.scatter_rev(std::plus<T>());
@@ -226,7 +226,7 @@ void solver(MPI_Comm comm)
   int num_it = linalg::cg(*u->x(), b, action, 200, 1e-6);
 
   // Set BC values in the solution vectors
-  fem::set_bc<T, U>(u->x()->mutable_array(), {bc}, T(1));
+  bc->set(u->x()->mutable_array(), std::nullopt, T(1));
 
   // Compute L2 error (squared) of the solution vector e = (u - u_d, u
   // - u_d)*dx
diff --git a/cpp/dolfinx/fem/DirichletBC.h b/cpp/dolfinx/fem/DirichletBC.h
index 89dd7d85677..b43066ff253 100644
--- a/cpp/dolfinx/fem/DirichletBC.h
+++ b/cpp/dolfinx/fem/DirichletBC.h
@@ -17,6 +17,7 @@
 #include <dolfinx/common/types.h>
 #include <functional>
 #include <memory>
+#include <optional>
 #include <span>
 #include <type_traits>
 #include <utility>
@@ -469,129 +470,118 @@ class DirichletBC
     return {_dofs0, _owned_indices0};
   }
 
-  /// Set bc entries in `x` to `scale * x_bc`
+  /// @brief Set entries in an array that are constrained by Dirichlet
+  /// boundary conditions.
   ///
-  /// @param[in] x The array in which to set `scale * x_bc[i]`, where
-  /// x_bc[i] is the boundary value of x[i]. Entries in x that do not
-  /// have a Dirichlet condition applied to them are unchanged. The
-  /// length of x must be less than or equal to the index of the
-  /// greatest boundary dof index. To set values only for
-  /// degrees-of-freedom that are owned by the calling rank, the length
-  /// of the array @p x should be equal to the number of dofs owned by
-  /// this rank.
-  /// @param[in] scale The scaling value to apply
-  void set(std::span<T> x, T scale = 1) const
+  /// Entries in `x` that are constrained by a Dirichlet boundary
+  /// conditions are set to `alpha * (x_bc - x0)`, where `x_bc` is the
+  /// (interpolated) boundary condition value.
+  ///
+  /// For elements with point-wise evaluated degrees-of-freedom, e.g.
+  /// Lagrange elements, `x_bc` is the value of the boundary condition
+  /// at the degree-of-freedom. For elements with moment
+  /// degrees-of-freedom, `x_bc` is the value of the boundary condition
+  /// interpolated into the finite element space.
+  ///
+  /// If `x` includes ghosted entries (entries available on the calling
+  /// rank but owned by another rank), ghosted entries constrained by a
+  /// Dirichlet condition will also be set.
+  ///
+  /// @param[in,out] x Array to modify for Dirichlet boundary
+  /// conditions.
+  /// @param[in] x0 Optional array used in computing the value to set.
+  /// If not provided it is treated as zero.
+  /// @param[in] alpha Scaling to apply.
+  void set(std::span<T> x, std::optional<std::span<const T>> x0,
+           T alpha = 1) const
   {
-    if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
+    std::int32_t x_size = x.size();
+    if (alpha == T(0)) // Optimisation for when alpha == 0
     {
-      auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
-      assert(g);
-      std::span<const T> values = g->x()->array();
-      auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
-      std::int32_t x_size = x.size();
-      for (std::size_t i = 0; i < _dofs0.size(); ++i)
+      for (std::int32_t idx : _dofs0)
       {
-        if (_dofs0[i] < x_size)
-        {
-          assert(dofs1_g[i] < (std::int32_t)values.size());
-          x[_dofs0[i]] = scale * values[dofs1_g[i]];
-        }
+        if (idx < x_size)
+          x[idx] = 0;
       }
     }
-    else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
-    {
-      auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
-      std::vector<T> value = g->value;
-      int bs = _function_space->dofmap()->bs();
-      std::int32_t x_size = x.size();
-      std::ranges::for_each(_dofs0,
-                            [x_size, bs, scale, &value, &x](auto dof)
-                            {
-                              if (dof < x_size)
-                                x[dof] = scale * value[dof % bs];
-                            });
-    }
-  }
-
-  /// Set bc entries in `x` to `scale * (x0 - x_bc)`
-  /// @param[in] x The array in which to set `scale * (x0 - x_bc)`
-  /// @param[in] x0 The array used in compute the value to set
-  /// @param[in] scale The scaling value to apply
-  void set(std::span<T> x, std::span<const T> x0, T scale = 1) const
-  {
-    if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
+    else
     {
-      auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
-      assert(g);
-      std::span<const T> values = g->x()->array();
-      assert(x.size() <= x0.size());
-      auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
-      std::int32_t x_size = x.size();
-      for (std::size_t i = 0; i < _dofs0.size(); ++i)
+      if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
       {
-        if (_dofs0[i] < x_size)
+        auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
+        assert(g);
+        auto dofs1_g
+            = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
+        std::span<const T> values = g->x()->array();
+        if (x0.has_value())
         {
-          assert(dofs1_g[i] < (std::int32_t)values.size());
-          x[_dofs0[i]] = scale * (values[dofs1_g[i]] - x0[_dofs0[i]]);
+          std::span<const T> _x0 = x0.value();
+          assert(x.size() <= _x0.size());
+          for (std::size_t i = 0; i < _dofs0.size(); ++i)
+          {
+            if (_dofs0[i] < x_size)
+            {
+              assert(dofs1_g[i] < (std::int32_t)values.size());
+              x[_dofs0[i]] = alpha * (values[dofs1_g[i]] - _x0[_dofs0[i]]);
+            }
+          }
+        }
+        else
+        {
+          for (std::size_t i = 0; i < _dofs0.size(); ++i)
+          {
+            if (_dofs0[i] < x_size)
+            {
+              assert(dofs1_g[i] < (std::int32_t)values.size());
+              x[_dofs0[i]] = alpha * values[dofs1_g[i]];
+            }
+          }
+        }
+      }
+      else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
+      {
+        auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
+        const std::vector<T>& value = g->value;
+        std::int32_t bs = _function_space->dofmap()->bs();
+        if (x0.has_value())
+        {
+          assert(x.size() <= x0.value().size());
+          std::ranges::for_each(
+              _dofs0,
+              [x_size, &x, x0 = x0.value(), &value, alpha, bs](auto dof)
+              {
+                if (dof < x_size)
+                  x[dof] = alpha * (value[dof % bs] - x0[dof]);
+              });
+        }
+        else
+        {
+          std::ranges::for_each(_dofs0,
+                                [x_size, bs, alpha, &value, &x](auto dof)
+                                {
+                                  if (dof < x_size)
+                                    x[dof] = alpha * value[dof % bs];
+                                });
         }
       }
-    }
-    else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
-    {
-      auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
-      const std::vector<T>& value = g->value;
-      std::int32_t bs = _function_space->dofmap()->bs();
-      std::ranges::for_each(_dofs0,
-                            [&x, &x0, &value, scale, bs](auto dof)
-                            {
-                              if (dof < (std::int32_t)x.size())
-                                x[dof] = scale * (value[dof % bs] - x0[dof]);
-                            });
     }
   }
 
-  /// @todo Review this function - it is almost identical to the
-  /// 'DirichletBC::set' function
+  /// @brief Set `markers[i] = true` if dof `i` has a boundary condition
+  /// applied.
   ///
-  /// Set boundary condition value for entries with an applied boundary
-  /// condition. Other entries are not modified.
-  /// @param[out] values The array in which to set the dof values.
-  /// The array must be at least as long as the array associated with V1
-  /// (the space of the function that provides the dof values)
-  void dof_values(std::span<T> values) const
-  {
-    if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
-    {
-      auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
-      assert(g);
-      std::span<const T> g_values = g->x()->array();
-      auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
-      for (std::size_t i = 0; i < dofs1_g.size(); ++i)
-        values[_dofs0[i]] = g_values[dofs1_g[i]];
-    }
-    else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
-    {
-      auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
-      assert(g);
-      const std::vector<T>& g_value = g->value;
-      const std::int32_t bs = _function_space->dofmap()->bs();
-      for (std::size_t i = 0; i < _dofs0.size(); ++i)
-        values[_dofs0[i]] = g_value[_dofs0[i] % bs];
-    }
-  }
-
-  /// Set markers[i] = true if dof i has a boundary condition applied.
-  /// Value of markers[i] is not changed otherwise.
-  /// @param[in,out] markers Entry makers[i] is set to true if dof i in
-  /// V0 had a boundary condition applied, i.e. dofs which are fixed by
-  /// a boundary condition. Other entries in @p markers are left
+  /// Value of `markers[i]` is not changed otherwise.
+  ///
+  /// @param[in,out] markers Entry `makers[i]` is set to true if dof `i`
+  /// in V0 had a boundary condition applied, i.e. dofs which are fixed
+  /// by a boundary condition. Other entries in `markers` are left
   /// unchanged.
   void mark_dofs(std::span<std::int8_t> markers) const
   {
-    for (std::size_t i = 0; i < _dofs0.size(); ++i)
+    for (std::int32_t idx : _dofs0)
     {
-      assert(_dofs0[i] < (std::int32_t)markers.size());
-      markers[_dofs0[i]] = true;
+      assert(idx < (std::int32_t)markers.size());
+      markers[idx] = true;
     }
   }
 
@@ -604,8 +594,8 @@ class DirichletBC
                std::shared_ptr<const Constant<T>>>
       _g;
 
-  // Dof indices (_dofs0) in _function_space and (_dofs1_g) in the
-  // space of _g. _dofs1_g may be empty if _dofs0 can be re-used
+  // Dof indices (_dofs0) in _function_space and (_dofs1_g) in the space
+  // of _g. _dofs1_g may be empty if _dofs0 can be re-used
   std::vector<std::int32_t> _dofs0, _dofs1_g;
 
   // The first _owned_indices in _dofs are owned by this process
diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h
index 6687a1acff4..7e5531de424 100644
--- a/cpp/dolfinx/fem/assemble_vector_impl.h
+++ b/cpp/dolfinx/fem/assemble_vector_impl.h
@@ -22,6 +22,7 @@
 #include <dolfinx/mesh/Topology.h>
 #include <functional>
 #include <memory>
+#include <optional>
 #include <span>
 #include <vector>
 
@@ -70,7 +71,7 @@ using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
 /// @param[in] bc_markers1 Marker to identify which DOFs have boundary
 /// conditions applied.
 /// @param[in] x0 Vector used in the lifting.
-/// @param[in] scale Scaling to apply.
+/// @param[in] alpha Scaling to apply.
 template <dolfinx::scalar T, int _bs0 = -1, int _bs1 = -1>
 void _lift_bc_cells(
     std::span<T> b, mdspan2_t x_dofmap,
@@ -83,7 +84,7 @@ void _lift_bc_cells(
     std::span<const T> coeffs, int cstride,
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T scale)
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
 {
   if (cells.empty())
     return;
@@ -183,9 +184,9 @@ void _lift_bc_cells(
             const T bc = bc_values1[jj];
             const T _x0 = x0.empty() ? 0 : x0[jj];
             // const T _x0 = 0;
-            // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
+            // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
             for (int m = 0; m < num_rows; ++m)
-              be[m] -= Ae[m * num_cols + _bs1 * j + k] * scale * (bc - _x0);
+              be[m] -= Ae[m * num_cols + _bs1 * j + k] * alpha * (bc - _x0);
           }
         }
       }
@@ -199,9 +200,9 @@ void _lift_bc_cells(
           {
             const T bc = bc_values1[jj];
             const T _x0 = x0.empty() ? 0 : x0[jj];
-            // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
+            // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
             for (int m = 0; m < num_rows; ++m)
-              be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
+              be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
           }
         }
       }
@@ -255,7 +256,7 @@ void _lift_bc_cells(
 /// @param[in] bc_markers1 Marker to identify which DOFs have boundary
 /// conditions applied.
 /// @param[in] x0 The vector used in the lifting.
-/// @param[in] scale The scaling to apply.
+/// @param[in] alpha The scaling to apply.
 /// @param[in] perms Facet permutation integer. Empty if facet
 /// permutations are not required.
 template <dolfinx::scalar T, int _bs = -1>
@@ -270,7 +271,7 @@ void _lift_bc_exterior_facets(
     std::span<const T> coeffs, int cstride,
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T scale,
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha,
     std::span<const std::uint8_t> perms)
 {
   if (facets.empty())
@@ -358,9 +359,9 @@ void _lift_bc_exterior_facets(
         {
           const T bc = bc_values1[jj];
           const T _x0 = x0.empty() ? 0 : x0[jj];
-          // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
+          // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
           for (int m = 0; m < num_rows; ++m)
-            be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
+            be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
         }
       }
     }
@@ -406,7 +407,7 @@ void _lift_bc_exterior_facets(
 /// @param[in] bc_markers1 Marker to identify which DOFs have boundary
 /// conditions applied.
 /// @param[in] x0 The vector used in the lifting.
-/// @param[in] scale The scaling to apply
+/// @param[in] alpha The scaling to apply
 template <dolfinx::scalar T, int _bs = -1>
 void _lift_bc_interior_facets(
     std::span<T> b, mdspan2_t x_dofmap,
@@ -420,7 +421,7 @@ void _lift_bc_interior_facets(
     std::span<const std::uint32_t> cell_info0,
     std::span<const std::uint32_t> cell_info1,
     std::span<const std::uint8_t> perms, std::span<const T> bc_values1,
-    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T scale)
+    std::span<const std::int8_t> bc_markers1, std::span<const T> x0, T alpha)
 {
   if (facets.empty())
     return;
@@ -568,9 +569,9 @@ void _lift_bc_interior_facets(
         {
           const T bc = bc_values1[jj];
           const T _x0 = x0.empty() ? 0 : x0[jj];
-          // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
+          // be -= Ae.col(bs1 * j + k) * alpha * (bc - _x0);
           for (int m = 0; m < num_rows; ++m)
-            be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
+            be[m] -= Ae[m * num_cols + bs1 * j + k] * alpha * (bc - _x0);
         }
       }
     }
@@ -586,11 +587,11 @@ void _lift_bc_interior_facets(
         {
           const T bc = bc_values1[jj];
           const T _x0 = x0.empty() ? 0 : x0[jj];
-          // be -= Ae.col(offset + bs1 * j + k) * scale * (bc - x0[jj]);
+          // be -= Ae.col(offset + bs1 * j + k) * alpha * (bc - x0[jj]);
           for (int m = 0; m < num_rows; ++m)
           {
             be[m]
-                -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0);
+                -= Ae[m * num_cols + offset + bs1 * j + k] * alpha * (bc - _x0);
           }
         }
       }
@@ -910,7 +911,7 @@ void assemble_interior_facets(
 
 /// Modify RHS vector to account for boundary condition such that:
 ///
-/// b <- b - scale * A (x_bc - x0)
+/// b <- b - alpha * A.(x_bc - x0)
 ///
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear form that generates A
@@ -923,7 +924,7 @@ void assemble_interior_facets(
 /// which bcs belong
 /// @param[in] x0 The array used in the lifting, typically a 'current
 /// solution' in a Newton method
-/// @param[in] scale Scaling to apply
+/// @param[in] alpha Scaling to apply
 template <dolfinx::scalar T, std::floating_point U>
 void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
              std::span<const scalar_value_type_t<T>> x,
@@ -932,7 +933,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
                             std::pair<std::span<const T>, int>>& coefficients,
              std::span<const T> bc_values1,
              std::span<const std::int8_t> bc_markers1, std::span<const T> x0,
-             T scale)
+             T alpha)
 {
   // Integration domain mesh
   std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
@@ -986,7 +987,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
           {dofmap0, bs0, a.domain(IntegralType::cell, i, *mesh0)}, P0,
           {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)}, P1T,
           constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
-          bc_markers1, x0, scale);
+          bc_markers1, x0, alpha);
     }
     else if (bs0 == 3 and bs1 == 3)
     {
@@ -995,7 +996,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
           {dofmap0, bs0, a.domain(IntegralType::cell, i, *mesh0)}, P0,
           {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)}, P1T,
           constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
-          bc_markers1, x0, scale);
+          bc_markers1, x0, alpha);
     }
     else
     {
@@ -1004,7 +1005,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
                      P0,
                      {dofmap1, bs1, a.domain(IntegralType::cell, i, *mesh1)},
                      P1T, constants, coeffs, cstride, cell_info0, cell_info1,
-                     bc_values1, bc_markers1, x0, scale);
+                     bc_values1, bc_markers1, x0, alpha);
     }
   }
 
@@ -1030,7 +1031,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
         {dofmap0, bs0, a.domain(IntegralType::exterior_facet, i, *mesh0)}, P0,
         {dofmap1, bs1, a.domain(IntegralType::exterior_facet, i, *mesh1)}, P1T,
         constants, coeffs, cstride, cell_info0, cell_info1, bc_values1,
-        bc_markers1, x0, scale, perms);
+        bc_markers1, x0, alpha, perms);
   }
 
   for (int i : a.integral_ids(IntegralType::interior_facet))
@@ -1045,19 +1046,20 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
         {dofmap0, bs0, a.domain(IntegralType::interior_facet, i, *mesh0)}, P0,
         {dofmap1, bs1, a.domain(IntegralType::interior_facet, i, *mesh1)}, P1T,
         constants, coeffs, cstride, cell_info0, cell_info1, perms, bc_values1,
-        bc_markers1, x0, scale);
+        bc_markers1, x0, alpha);
   }
 }
 
 /// Modify b such that:
 ///
-///   b <- b - scale * A_j (g_j - x0_j)
+///   b <- b - alpha * A_j.(g_j - x0_j)
+///
+/// where j is a block (nest) row index. For a non-blocked problem j =
+/// 0. The boundary conditions bc1 are on the trial spaces V_j. The
+/// forms in [a] must have the same test space as L (from which b was
+/// built), but the trial space may differ. If x0 is not supplied, then
+/// it is treated as zero.
 ///
-/// where j is a block (nest) row index. For a non-blocked problem j = 0.
-/// The boundary conditions bc1 are on the trial spaces V_j. The forms
-/// in [a] must have the same test space as L (from which b was built),
-/// but the trial space may differ. If x0 is not supplied, then it is
-/// treated as zero.
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear forms, where a[j] is the form that
 /// generates A_j
@@ -1067,7 +1069,7 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
 /// bcs1[2] are the boundary conditions applied to the columns of a[2] /
 /// x0[2] block
 /// @param[in] x0 The vectors used in the lifting
-/// @param[in] scale Scaling to apply
+/// @param[in] alpha Scaling to apply
 template <dolfinx::scalar T, std::floating_point U>
 void apply_lifting(
     std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>> a,
@@ -1076,7 +1078,7 @@ void apply_lifting(
                                std::pair<std::span<const T>, int>>>& coeffs,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
         bcs1,
-    const std::vector<std::span<const T>>& x0, T scale)
+    const std::vector<std::span<const T>>& x0, T alpha)
 {
   if (!x0.empty() and x0.size() != a.size())
   {
@@ -1116,31 +1118,32 @@ void apply_lifting(
       for (const std::shared_ptr<const DirichletBC<T, U>>& bc : bcs1[j])
       {
         bc->mark_dofs(bc_markers1);
-        bc->dof_values(bc_values1);
+        // bc->dof_values(bc_values1);
+        bc->set(bc_values1, std::nullopt, 1);
       }
 
       if (!x0.empty())
       {
         lift_bc<T>(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1,
-                   bc_markers1, x0[j], scale);
+                   bc_markers1, x0[j], alpha);
       }
       else
       {
         lift_bc<T>(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1,
-                   bc_markers1, std::span<const T>(), scale);
+                   bc_markers1, std::span<const T>(), alpha);
       }
     }
   }
 }
 
-/// Assemble linear form into a vector
+/// @brief Assemble linear form into a vector.
 /// @param[in,out] b The vector to be assembled. It will not be zeroed
 /// before assembly.
-/// @param[in] L The linear forms to assemble into b
-/// @param[in] x_dofmap Mesh geometry dofmap
-/// @param[in] x Mesh coordinates
-/// @param[in] constants Packed constants that appear in `L`
-/// @param[in] coefficients Packed coefficients that appear in `L`
+/// @param[in] L Linear forms to assemble into b.
+/// @param[in] x_dofmap Mesh geometry dofmap.
+/// @param[in] x Mesh coordinates.
+/// @param[in] constants Packed constants that appear in `L`.
+/// @param[in] coefficients Packed coefficients that appear in `L`.
 template <dolfinx::scalar T, std::floating_point U>
 void assemble_vector(
     std::span<T> b, const Form<T, U>& L, mdspan2_t x_dofmap,
diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
index c53aef84b40..decd5d35411 100644
--- a/cpp/dolfinx/fem/assembler.h
+++ b/cpp/dolfinx/fem/assembler.h
@@ -15,6 +15,7 @@
 #include <cstdint>
 #include <dolfinx/common/types.h>
 #include <memory>
+#include <optional>
 #include <span>
 #include <vector>
 
@@ -37,10 +38,10 @@ make_coefficients_span(const std::map<std::pair<IntegralType, int>,
 {
   using Key = typename std::remove_reference_t<decltype(coeffs)>::key_type;
   std::map<Key, std::pair<std::span<const T>, int>> c;
-  std::ranges::transform(coeffs, std::inserter(c, c.end()),
-                         [](auto& e) -> typename decltype(c)::value_type {
-                           return {e.first, {e.second.first, e.second.second}};
-                         });
+  std::ranges::transform(
+      coeffs, std::inserter(c, c.end()),
+      [](auto& e) -> typename decltype(c)::value_type
+      { return {e.first, {e.second.first, e.second.second}}; });
   return c;
 }
 
@@ -137,7 +138,7 @@ void assemble_vector(std::span<T> b, const Form<T, U>& L)
 
 /// Modify b such that:
 ///
-///   b <- b - scale * A_j (g_j - x0_j)
+///   b <- b - alpha * A_j (g_j - x0_j)
 ///
 /// where j is a block (nest) index. For a non-blocked problem j = 0. The
 /// boundary conditions bcs1 are on the trial spaces V_j. The forms in
@@ -155,18 +156,18 @@ void apply_lifting(
                                std::pair<std::span<const T>, int>>>& coeffs,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
         bcs1,
-    const std::vector<std::span<const T>>& x0, T scale)
+    const std::vector<std::span<const T>>& x0, T alpha)
 {
   // If all forms are null, there is nothing to do
   if (std::ranges::all_of(a, [](auto ptr) { return ptr == nullptr; }))
     return;
 
-  impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, scale);
+  impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, alpha);
 }
 
 /// Modify b such that:
 ///
-///   b <- b - scale * A_j (g_j - x0_j)
+///   b <- b - alpha * A_j.(g_j - x0_j)
 ///
 /// where j is a block (nest) index. For a non-blocked problem j = 0. The
 /// boundary conditions bcs1 are on the trial spaces V_j. The forms in
@@ -181,7 +182,7 @@ void apply_lifting(
     std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>>& a,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
         bcs1,
-    const std::vector<std::span<const T>>& x0, T scale)
+    const std::vector<std::span<const T>>& x0, T alpha)
 {
   std::vector<
       std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
@@ -211,7 +212,7 @@ void apply_lifting(
   std::ranges::transform(coeffs, std::back_inserter(_coeffs),
                          [](auto& c) { return make_coefficients_span(c); });
 
-  apply_lifting(b, a, _constants, _coeffs, bcs1, x0, scale);
+  apply_lifting(b, a, _constants, _coeffs, bcs1, x0, alpha);
 }
 
 // -- Matrices ---------------------------------------------------------------
@@ -404,44 +405,4 @@ void set_diagonal(
     }
   }
 }
-
-// -- Setting bcs ------------------------------------------------------------
-
-// FIXME: Move these function elsewhere?
-
-// FIXME: clarify x0
-// FIXME: clarify what happens with ghosts
-
-/// Set bc values in owned (local) part of the vector, multiplied by
-/// 'scale'. The vectors b and x0 must have the same local size. The bcs
-/// should be on (sub-)spaces of the form L that b represents.
-template <dolfinx::scalar T, std::floating_point U>
-void set_bc(std::span<T> b,
-            const std::vector<std::shared_ptr<const DirichletBC<T, U>>>& bcs,
-            std::span<const T> x0, T scale = 1)
-{
-  if (b.size() > x0.size())
-    throw std::runtime_error("Size mismatch between b and x0 vectors.");
-  for (auto& bc : bcs)
-  {
-    assert(bc);
-    bc->set(b, x0, scale);
-  }
-}
-
-/// Set bc values in owned (local) part of the vector, multiplied by
-/// 'scale'. The bcs should be on (sub-)spaces of the form L that b
-/// represents.
-template <dolfinx::scalar T, std::floating_point U>
-void set_bc(std::span<T> b,
-            const std::vector<std::shared_ptr<const DirichletBC<T, U>>>& bcs,
-            T scale = 1)
-{
-  for (auto& bc : bcs)
-  {
-    assert(bc);
-    bc->set(b, scale);
-  }
-}
-
 } // namespace dolfinx::fem
diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h
index edb4aa4bfdc..a94af0b0513 100644
--- a/cpp/dolfinx/fem/petsc.h
+++ b/cpp/dolfinx/fem/petsc.h
@@ -34,7 +34,7 @@ class DirichletBC;
 /// @brief Helper functions for assembly into PETSc data structures
 namespace petsc
 {
-/// Create a matrix
+/// @brief Create a matrix
 /// @param[in] a A bilinear form
 /// @param[in] type The PETSc matrix type to create
 /// @return A sparse matrix with a layout and sparsity that matches the
@@ -49,7 +49,8 @@ Mat create_matrix(const Form<PetscScalar, T>& a,
   return la::petsc::create_matrix(a.mesh()->comm(), pattern, type);
 }
 
-/// Initialise a monolithic matrix for an array of bilinear forms
+/// @brief Initialise a monolithic matrix for an array of bilinear
+/// forms.
 /// @param[in] a Rectangular array of bilinear forms. The `a(i, j)` form
 /// will correspond to the `(i, j)` block in the returned matrix
 /// @param[in] type The type of PETSc Mat. If empty the PETSc default is
@@ -316,7 +317,7 @@ void assemble_vector(Vec b, const Form<PetscScalar, T>& L)
 ///
 /// Modify b such that:
 ///
-///   b <- b - scale * A_j (g_j - x0_j)
+///   b <- b - alpha * A_j (g_j - x0_j)
 ///
 /// where j is a block (nest) index. For a non-blocked problem j = 0. The
 /// boundary conditions bcs1 are on the trial spaces V_j. The forms in
@@ -335,7 +336,7 @@ void apply_lifting(
         coeffs,
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar, T>>>>& bcs1,
-    const std::vector<Vec>& x0, PetscScalar scale)
+    const std::vector<Vec>& x0, PetscScalar alpha)
 {
   Vec b_local;
   VecGhostGetLocalForm(b, &b_local);
@@ -346,7 +347,7 @@ void apply_lifting(
   std::span<PetscScalar> _b(array, n);
 
   if (x0.empty())
-    fem::apply_lifting(_b, a, constants, coeffs, bcs1, {}, scale);
+    fem::apply_lifting(_b, a, constants, coeffs, bcs1, {}, alpha);
   else
   {
     std::vector<std::span<const PetscScalar>> x0_ref;
@@ -363,7 +364,7 @@ void apply_lifting(
     }
 
     std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
-    fem::apply_lifting(_b, a, constants, coeffs, bcs1, x0_tmp, scale);
+    fem::apply_lifting(_b, a, constants, coeffs, bcs1, x0_tmp, alpha);
 
     for (std::size_t i = 0; i < x0_local.size(); ++i)
     {
@@ -384,7 +385,7 @@ void apply_lifting(
 
 /// Modify b such that:
 ///
-///   b <- b - scale * A_j (g_j - x0_j)
+///   b <- b - alpha * A_j (g_j - x0_j)
 ///
 /// where j is a block (nest) index. For a non-blocked problem j = 0. The
 /// boundary conditions bcs1 are on the trial spaces V_j. The forms in
@@ -401,7 +402,7 @@ void apply_lifting(
     const std::vector<
         std::vector<std::shared_ptr<const DirichletBC<PetscScalar, double>>>>&
         bcs1,
-    const std::vector<Vec>& x0, PetscScalar scale)
+    const std::vector<Vec>& x0, PetscScalar alpha)
 {
   Vec b_local;
   VecGhostGetLocalForm(b, &b_local);
@@ -412,7 +413,7 @@ void apply_lifting(
   std::span<PetscScalar> _b(array, n);
 
   if (x0.empty())
-    fem::apply_lifting<PetscScalar>(_b, a, bcs1, {}, scale);
+    fem::apply_lifting<PetscScalar>(_b, a, bcs1, {}, alpha);
   else
   {
     std::vector<std::span<const PetscScalar>> x0_ref;
@@ -429,7 +430,7 @@ void apply_lifting(
     }
 
     std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
-    fem::apply_lifting<PetscScalar>(_b, a, bcs1, x0_tmp, scale);
+    fem::apply_lifting<PetscScalar>(_b, a, bcs1, x0_tmp, alpha);
 
     for (std::size_t i = 0; i < x0_local.size(); ++i)
     {
@@ -450,13 +451,13 @@ void apply_lifting(
 // FIXME: clarify what happens with ghosts
 
 /// Set bc values in owned (local) part of the PETSc vector, multiplied
-/// by 'scale'. The vectors b and x0 must have the same local size. The
+/// by 'alpha'. The vectors b and x0 must have the same local size. The
 /// bcs should be on (sub-)spaces of the form L that b represents.
 template <std::floating_point T>
 void set_bc(
     Vec b,
     const std::vector<std::shared_ptr<const DirichletBC<PetscScalar, T>>>& bcs,
-    const Vec x0, PetscScalar scale = 1)
+    const Vec x0, PetscScalar alpha = 1)
 {
   PetscInt n = 0;
   VecGetLocalSize(b, &n);
@@ -472,13 +473,16 @@ void set_bc(
     const PetscScalar* array = nullptr;
     VecGetArrayRead(x0_local, &array);
     std::span<const PetscScalar> _x0(array, n);
-    fem::set_bc(_b, bcs, _x0, scale);
+    for (auto& bc : bcs)
+      bc->set(_b, _x0, alpha);
     VecRestoreArrayRead(x0_local, &array);
     VecGhostRestoreLocalForm(x0, &x0_local);
   }
   else
-    fem::set_bc(_b, bcs, scale);
-
+  {
+    for (auto& bc : bcs)
+      bc->set(_b, std::nullopt, alpha);
+  }
   VecRestoreArray(b, &array);
 }
 
diff --git a/python/demo/demo_elasticity.py b/python/demo/demo_elasticity.py
index 79c4856a6ff..a9f69887a17 100644
--- a/python/demo/demo_elasticity.py
+++ b/python/demo/demo_elasticity.py
@@ -50,7 +50,7 @@
     functionspace,
     locate_dofs_topological,
 )
-from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
+from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import CellType, GhostMode, create_box, locate_entities_boundary
 from ufl import dx, grad, inner
@@ -183,7 +183,7 @@ def σ(v):
 b = assemble_vector(L)
 apply_lifting(b, [a], bcs=[[bc]])
 b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)  # type: ignore
-set_bc(b, [bc])
+bc.set(b.array_w)
 # -
 
 # Create the near-nullspace and attach it to the PETSc matrix:
diff --git a/python/demo/demo_poisson_matrix_free.py b/python/demo/demo_poisson_matrix_free.py
index eaeb417d310..3d5eb8fc638 100644
--- a/python/demo/demo_poisson_matrix_free.py
+++ b/python/demo/demo_poisson_matrix_free.py
@@ -164,12 +164,12 @@
 # Apply lifting: b <- b - A * x_bc
 b = fem.assemble_vector(L_fem)
 ui.x.array[:] = 0.0
-fem.set_bc(ui.x.array, [bc], scale=-1.0)
+bc.set(ui.x.array, alpha=-1.0)
 fem.assemble_vector(b.array, M_fem)
 b.scatter_reverse(la.InsertMode.add)
 
 # Set BC dofs to zero on right hand side
-fem.set_bc(b.array, [bc], scale=0.0)
+bc.set(b.array, alpha=0.0)
 b.scatter_forward()
 
 # To implement the matrix-free CG solver using *DOLFINx* vectors, we define the
@@ -188,7 +188,7 @@ def action_A(x, y):
     y.scatter_reverse(la.InsertMode.add)
 
     # Set BC dofs to zero
-    fem.set_bc(y.array, [bc], scale=0.0)
+    bc.set(y.array, alpha=0.0)
 
 
 # ### Basic conjugate gradient solver
@@ -248,7 +248,7 @@ def _global_dot(comm, v0, v1):
 iter_cg1 = cg(mesh.comm, action_A, u.x, b, max_iter=200, rtol=rtol)
 
 # Set BC values in the solution vector
-fem.set_bc(u.x.array, [bc], scale=1.0)
+bc.set(u.x.array, alpha=1.0)
 
 
 def L2Norm(u):
diff --git a/python/demo/demo_pyamg.py b/python/demo/demo_pyamg.py
index f59f8588c6d..69de02486fe 100644
--- a/python/demo/demo_pyamg.py
+++ b/python/demo/demo_pyamg.py
@@ -46,7 +46,6 @@
     form,
     functionspace,
     locate_dofs_topological,
-    set_bc,
 )
 from dolfinx.mesh import CellType, create_box, locate_entities_boundary
 from ufl import ds, dx, grad, inner
@@ -95,7 +94,7 @@ def poisson_problem(dtype: npt.DTypeLike, solver_type: str):
     A = assemble_matrix(a, [bc]).to_scipy()
     b = assemble_vector(L)
     apply_lifting(b.array, [a], bcs=[[bc]])
-    set_bc(b.array, [bc])
+    bc.set(b.array)
 
     # Create solution function and create AMG solver
     uh = fem.Function(V, dtype=dtype)
@@ -208,7 +207,7 @@ def σ(v):
     A = assemble_matrix(a, bcs=[bc]).to_scipy()
     b = assemble_vector(L)
     apply_lifting(b.array, [a], bcs=[[bc]])
-    set_bc(b.array, [bc])
+    bc.set(b.array)
 
     uh = fem.Function(V, dtype=dtype)
     B = nullspace_elasticty(V)
diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py
index 93c7d29d690..953f7f64af2 100644
--- a/python/demo/demo_static-condensation.py
+++ b/python/demo/demo_static-condensation.py
@@ -55,7 +55,7 @@
     functionspace,
     locate_dofs_topological,
 )
-from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
+from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector
 from dolfinx.io import XDMFFile
 from dolfinx.jit import ffcx_jit
 from dolfinx.mesh import locate_entities_boundary, meshtags
@@ -170,7 +170,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL):
 b = assemble_vector(b1)
 apply_lifting(b, [a_cond], bcs=[[bc]])
 b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)  # type: ignore
-set_bc(b, [bc])
+bc.set(b)
 
 uc = Function(U)
 solver = PETSc.KSP().create(A_cond.getComm())  # type: ignore
diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py
index 65de072c569..ab14e9b3437 100644
--- a/python/demo/demo_stokes.py
+++ b/python/demo/demo_stokes.py
@@ -521,7 +521,8 @@ def mixed_direct():
     b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
 
     # Set Dirichlet boundary condition values in the RHS
-    fem.petsc.set_bc(b, bcs)
+    for bc in bcs:
+        bc.set(b)
 
     # Create and configure solver
     ksp = PETSc.KSP().create(msh.comm)
diff --git a/python/demo/demo_types.py b/python/demo/demo_types.py
index f5914368a79..7d647ecfde9 100644
--- a/python/demo/demo_types.py
+++ b/python/demo/demo_types.py
@@ -141,7 +141,7 @@ def poisson(dtype):
     b = fem.assemble_vector(L0)
     fem.apply_lifting(b.array, [a0], bcs=[[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b.array, [bc])
+    bc.set(b.array)
 
     # Create a SciPy CSR  matrix that shares data with A and solve
     As = A.to_scipy()
@@ -206,7 +206,7 @@ def σ(v):
     b = fem.assemble_vector(L0)
     fem.apply_lifting(b.array, [a0], bcs=[[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b.array, [bc])
+    bc.set(b.array)
 
     # Create a SciPy CSR  matrix that shares data with A and solve
     As = A.to_scipy()
diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py
index 52a865219cd..be186a29050 100644
--- a/python/dolfinx/fem/__init__.py
+++ b/python/dolfinx/fem/__init__.py
@@ -173,7 +173,6 @@ def compute_integration_domains(
     "assemble_matrix",
     "assemble_vector",
     "apply_lifting",
-    "set_bc",
     "DirichletBC",
     "dirichletbc",
     "bcs_by_block",
@@ -192,4 +191,5 @@ def compute_integration_domains(
     "form_cpp_class",
     "create_form",
     "compile_form",
+    "set_bc",
 ]
diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py
index d3bb918e808..7b8f6eae6d3 100644
--- a/python/dolfinx/fem/assemble.py
+++ b/python/dolfinx/fem/assemble.py
@@ -10,6 +10,7 @@
 import collections
 import functools
 import typing
+import warnings
 
 import numpy as np
 
@@ -305,7 +306,7 @@ def apply_lifting(
     a: list[Form],
     bcs: list[list[DirichletBC]],
     x0: typing.Optional[list[np.ndarray]] = None,
-    scale: float = 1.0,
+    alpha: float = 1,
     constants=None,
     coeffs=None,
 ) -> None:
@@ -344,23 +345,27 @@ def apply_lifting(
     )
     _a = [None if form is None else form._cpp_object for form in a]
     _bcs = [[bc._cpp_object for bc in bcs0] for bcs0 in bcs]
-    _cpp.fem.apply_lifting(b, _a, constants, coeffs, _bcs, x0, scale)
+    _cpp.fem.apply_lifting(b, _a, constants, coeffs, _bcs, x0, alpha)
 
 
 def set_bc(
     b: np.ndarray,
     bcs: list[DirichletBC],
     x0: typing.Optional[np.ndarray] = None,
-    scale: float = 1.0,
+    scale: float = 1,
 ) -> None:
     """Insert boundary condition values into vector.
 
+    Note:
+        This function is deprecated.
+
     Only local (owned) entries are set, hence communication after
     calling this function is not required unless ghost entries need to
     be updated to the boundary condition value.
     """
-    _bcs = [bc._cpp_object for bc in bcs]
-    if x0 is None:
-        _cpp.fem.set_bc(b, _bcs, scale)
-    else:
-        _cpp.fem.set_bc(b, _bcs, x0, scale)
+    warnings.warn(
+        "dolfinx.fem.assembler.set_bc is deprecated. Use dolfinx.fem.DirichletBC.set instead.",
+        DeprecationWarning,
+    )
+    for bc in bcs:
+        bc.set(b, x0, scale)
diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py
index 737ca62d1e4..e433bc81d95 100644
--- a/python/dolfinx/fem/bcs.py
+++ b/python/dolfinx/fem/bcs.py
@@ -11,7 +11,7 @@
 import numbers
 import typing
 
-import numpy.typing
+import numpy.typing as npt
 
 if typing.TYPE_CHECKING:
     from dolfinx.fem.function import Constant, Function
@@ -57,7 +57,7 @@ def locate_dofs_geometrical(
 def locate_dofs_topological(
     V: typing.Union[dolfinx.fem.FunctionSpace, typing.Iterable[dolfinx.fem.FunctionSpace]],
     entity_dim: int,
-    entities: numpy.typing.NDArray[np.int32],
+    entities: npt.NDArray[np.int32],
     remote: bool = True,
 ) -> np.ndarray:
     """Locate degrees-of-freedom belonging to mesh entities topologically.
@@ -132,7 +132,34 @@ def function_space(self) -> dolfinx.fem.FunctionSpace:
         """The function space on which the boundary condition is defined"""
         return self._cpp_object.function_space
 
-    def dof_indices(self) -> tuple[np.ndarray, int]:
+    def set(
+        self, x: npt.NDArray, x0: typing.Optional[npt.NDArray[np.int32]] = None, alpha: float = 1
+    ) -> None:
+        """Set entries in an array that are constrained by Dirichlet boundary conditions.
+
+        Entries in ``x`` that are constrained by a Dirichlet boundary
+        conditions are set to ``alpha * (x_bc - x0)``, where ``x_bc`` is
+        the (interpolated) boundary condition value.
+
+        For elements with point-wise evaluated degrees-of-freedom, e.g.
+        Lagrange elements, ``x_bc`` is the value of the boundary condition
+        at the degree-of-freedom. For elements with moment
+        degrees-of-freedom, ``x_bc`` is the value of the boundary condition
+        interpolated into the finite element space.
+
+        If `x` includes ghosted entries (entries available on the calling
+        rank but owned by another rank), ghosted entries constrained by a
+        Dirichlet condition will also be set.
+
+        Args:
+            x: Array to modify for Dirichlet boundary conditions.
+            x0: Optional array used in computing the value to set. If
+                not provided it is treated as zero.
+            alpha: Scaling factor.
+        """
+        self._cpp_object.set(x, x0, alpha)
+
+    def dof_indices(self) -> tuple[npt.NDArray[np.int32], int]:
         """Access dof indices `(local indices, unrolled)`, including ghosts, to
         which a Dirichlet condition is applied, and the index to the first
         non-owned (ghost) index. The array of indices is sorted.
@@ -150,7 +177,7 @@ def dof_indices(self) -> tuple[np.ndarray, int]:
 
 def dirichletbc(
     value: typing.Union[Function, Constant, np.ndarray],
-    dofs: numpy.typing.NDArray[np.int32],
+    dofs: npt.NDArray[np.int32],
     V: typing.Optional[dolfinx.fem.FunctionSpace] = None,
 ) -> DirichletBC:
     """Create a representation of Dirichlet boundary condition which
diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py
index b71f0416ce2..02e3d84a469 100644
--- a/python/dolfinx/fem/petsc.py
+++ b/python/dolfinx/fem/petsc.py
@@ -305,7 +305,7 @@ def assemble_vector_block(
     a: list[list[Form]],
     bcs: list[DirichletBC] = [],
     x0: typing.Optional[PETSc.Vec] = None,
-    scale: float = 1.0,
+    alpha: float = 1,
     constants_L=None,
     coeffs_L=None,
     constants_a=None,
@@ -323,7 +323,7 @@ def assemble_vector_block(
     with b.localForm() as b_local:
         b_local.set(0.0)
     return _assemble_vector_block_vec(
-        b, L, a, bcs, x0, scale, constants_L, coeffs_L, constants_a, coeffs_a
+        b, L, a, bcs, x0, alpha, constants_L, coeffs_L, constants_a, coeffs_a
     )
 
 
@@ -334,7 +334,7 @@ def _assemble_vector_block_vec(
     a: list[list[Form]],
     bcs: list[DirichletBC] = [],
     x0: typing.Optional[PETSc.Vec] = None,
-    scale: float = 1.0,
+    alpha: float = 1,
     constants_L=None,
     coeffs_L=None,
     constants_a=None,
@@ -398,7 +398,7 @@ def _assemble_vector_block_vec(
     ):
         _cpp.fem.assemble_vector(b_sub, L_sub._cpp_object, const_L, coeff_L)
         _a_sub = [None if form is None else form._cpp_object for form in a_sub]
-        _cpp.fem.apply_lifting(b_sub, _a_sub, const_a, coeff_a, bcs1, x0_local, scale)
+        _cpp.fem.apply_lifting(b_sub, _a_sub, const_a, coeff_a, bcs1, x0_local, alpha)
 
     _cpp.la.petsc.scatter_local_vectors(b, b_local, maps)
     b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
@@ -406,12 +406,10 @@ def _assemble_vector_block_vec(
     bcs0 = _bcs_by_block(_extract_spaces(L), _bcs)
     offset = 0
     b_array = b.getArray(readonly=False)
-    for submap, bc, _x0 in zip(maps, bcs0, x0_sub):
+    for submap, bcs, _x0 in zip(maps, bcs0, x0_sub):
         size = submap[0].size_local * submap[1]
-        if _x0 is None:
-            _cpp.fem.set_bc(b_array[offset : offset + size], bc, scale)
-        else:
-            _cpp.fem.set_bc(b_array[offset : offset + size], bc, _x0, scale)
+        for bc in bcs:
+            bc.set(b_array[offset : offset + size], _x0, alpha)
         offset += size
 
     return b
@@ -434,8 +432,9 @@ def assemble_matrix(
     Args:
         a: Bilinear form to assembled into a matrix.
         bc: Dirichlet boundary conditions applied to the system.
-        diagonal: Value to set on the matrix diagonal for Dirichlet boundary
-            condition constrained degrees-of-freedom belonging to the same trial and test space.
+        diagonal: Value to set on the matrix diagonal for Dirichlet
+            boundary condition constrained degrees-of-freedom belonging
+            to the same trial and test space.
         constants: Constants appearing the in the form.
         coeffs: Coefficients appearing the in the form.
 
@@ -488,8 +487,9 @@ def assemble_matrix_nest(
         a: Rectangular (list-of-lists) array for bilinear forms.
         bcs: Dirichlet boundary conditions.
         mat_types: PETSc matrix type for each matrix block.
-        diagonal: Value to set on the matrix diagonal for Dirichlet boundary
-            condition constrained degrees-of-freedom belonging to the same trial and test space.
+        diagonal: Value to set on the matrix diagonal for Dirichlet
+            boundary condition constrained degrees-of-freedom belonging
+            to the same trial and test space.
         constants: Constants appearing the in the form.
         coeffs: Coefficients appearing the in the form.
 
@@ -520,8 +520,9 @@ def _assemble_matrix_nest_mat(
         a: Rectangular (list-of-lists) array for bilinear forms.
         bcs: Dirichlet boundary conditions.
         mat_types: PETSc matrix type for each matrix block.
-        diagonal: Value to set on the matrix diagonal for Dirichlet boundary
-            condition constrained degrees-of-freedom belonging to the same trial and test space.
+        diagonal: Value to set on the matrix diagonal for Dirichlet
+            boundary condition constrained degrees-of-freedom belonging
+            to the same trial and test space.
         constants: Constants appearing the in the form.
         coeffs: Coefficients appearing the in the form.
 
@@ -659,7 +660,7 @@ def apply_lifting(
     a: list[Form],
     bcs: list[list[DirichletBC]],
     x0: list[PETSc.Vec] = [],
-    scale: float = 1.0,
+    alpha: float = 1,
     constants=None,
     coeffs=None,
 ) -> None:
@@ -668,7 +669,7 @@ def apply_lifting(
         x0 = [stack.enter_context(x.localForm()) for x in x0]
         x0_r = [x.array_r for x in x0]
         b_local = stack.enter_context(b.localForm())
-        _assemble.apply_lifting(b_local.array_w, a, bcs, x0_r, scale, constants, coeffs)
+        _assemble.apply_lifting(b_local.array_w, a, bcs, x0_r, alpha, constants, coeffs)
 
 
 def apply_lifting_nest(
@@ -676,13 +677,12 @@ def apply_lifting_nest(
     a: list[list[Form]],
     bcs: list[DirichletBC],
     x0: typing.Optional[PETSc.Vec] = None,
-    scale: float = 1.0,
+    alpha: float = 1,
     constants=None,
     coeffs=None,
 ) -> PETSc.Vec:
     """Apply the function :func:`dolfinx.fem.apply_lifting` to each sub-vector
-    in a nested PETSc Vector.
-    """
+    in a nested PETSc Vector."""
     x0 = [] if x0 is None else x0.getNestSubVecs()
     bcs1 = _bcs_by_block(_extract_spaces(a, 1), bcs)
     constants = (
@@ -707,30 +707,31 @@ def apply_lifting_nest(
         else coeffs
     )
     for b_sub, a_sub, const, coeff in zip(b.getNestSubVecs(), a, constants, coeffs):
-        apply_lifting(b_sub, a_sub, bcs1, x0, scale, const, coeff)
+        apply_lifting(b_sub, a_sub, bcs1, x0, alpha, const, coeff)
     return b
 
 
 def set_bc(
-    b: PETSc.Vec, bcs: list[DirichletBC], x0: typing.Optional[PETSc.Vec] = None, scale: float = 1.0
+    b: PETSc.Vec, bcs: list[DirichletBC], x0: typing.Optional[PETSc.Vec] = None, alpha: float = 1
 ) -> None:
     """Apply the function :func:`dolfinx.fem.set_bc` to a PETSc Vector."""
     if x0 is not None:
         x0 = x0.array_r
-    _assemble.set_bc(b.array_w, bcs, x0, scale)
+    for bc in bcs:
+        bc.set(b.array_w, x0, alpha)
 
 
 def set_bc_nest(
     b: PETSc.Vec,
     bcs: list[list[DirichletBC]],
     x0: typing.Optional[PETSc.Vec] = None,
-    scale: float = 1.0,
+    alpha: float = 1,
 ) -> None:
     """Apply the function :func:`dolfinx.fem.set_bc` to each sub-vector of a nested PETSc Vector."""
     _b = b.getNestSubVecs()
     x0 = len(_b) * [None] if x0 is None else x0.getNestSubVecs()
     for b_sub, bc, x_sub in zip(_b, bcs, x0):
-        set_bc(b_sub, bc, x_sub, scale)
+        set_bc(b_sub, bc, x_sub, alpha)
 
 
 class LinearProblem:
@@ -848,7 +849,8 @@ def solve(self) -> _Function:
         # Apply boundary conditions to the rhs
         apply_lifting(self._b, [self._a], bcs=[self.bcs])
         self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
-        set_bc(self._b, self.bcs)
+        for bc in self.bcs:
+            bc.set(self._b.array_w)
 
         # Solve linear system and update ghost values in the solution
         self._solver.solve(self._b, self._x)
@@ -964,7 +966,7 @@ def F(self, x: PETSc.Vec, b: PETSc.Vec) -> None:
         assemble_vector(b, self._L)
 
         # Apply boundary condition
-        apply_lifting(b, [self._a], bcs=[self.bcs], x0=[x], scale=-1.0)
+        apply_lifting(b, [self._a], bcs=[self.bcs], x0=[x], alpha=-1.0)
         b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
         set_bc(b, self.bcs, x, -1.0)
 
diff --git a/python/dolfinx/wrappers/assemble.cpp b/python/dolfinx/wrappers/assemble.cpp
index d216039d083..6e9c04689a1 100644
--- a/python/dolfinx/wrappers/assemble.cpp
+++ b/python/dolfinx/wrappers/assemble.cpp
@@ -30,6 +30,7 @@
 #include <nanobind/stl/complex.h>
 #include <nanobind/stl/function.h>
 #include <nanobind/stl/map.h>
+#include <nanobind/stl/optional.h>
 #include <nanobind/stl/pair.h>
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/tuple.h>
@@ -193,7 +194,8 @@ void declare_assembly_functions(nb::module_& m)
       nb::arg("form"), "Pack coefficients for a Form.");
   m.def(
       "pack_constants",
-      [](const dolfinx::fem::Form<T, U>& form) {
+      [](const dolfinx::fem::Form<T, U>& form)
+      {
         return dolfinx_wrappers::as_nbarray(dolfinx::fem::pack_constants(form));
       },
       nb::arg("form"), "Pack constants for a Form.");
@@ -389,7 +391,7 @@ void declare_assembly_functions(nb::module_& m)
          const std::vector<std::vector<
              std::shared_ptr<const dolfinx::fem::DirichletBC<T, U>>>>& bcs1,
          const std::vector<nb::ndarray<const T, nb::ndim<1>, nb::c_contig>>& x0,
-         T scale)
+         T alpha)
       {
         std::vector<std::span<const T>> _x0;
         for (auto x : x0)
@@ -407,31 +409,11 @@ void declare_assembly_functions(nb::module_& m)
                                [](auto& c) { return py_to_cpp_coeffs(c); });
 
         dolfinx::fem::apply_lifting<T>(std::span<T>(b.data(), b.size()), a,
-                                       _constants, _coeffs, bcs1, _x0, scale);
+                                       _constants, _coeffs, bcs1, _x0, alpha);
       },
       nb::arg("b").noconvert(), nb::arg("a"), nb::arg("constants"),
-      nb::arg("coeffs"), nb::arg("bcs1"), nb::arg("x0"), nb::arg("scale"),
+      nb::arg("coeffs"), nb::arg("bcs1"), nb::arg("x0"), nb::arg("alpha"),
       "Modify vector for lifted boundary conditions");
-  m.def(
-      "set_bc",
-      [](nb::ndarray<T, nb::ndim<1>, nb::c_contig> b,
-         const std::vector<
-             std::shared_ptr<const dolfinx::fem::DirichletBC<T, U>>>& bcs,
-         nb::ndarray<const T, nb::ndim<1>, nb::c_contig> x0, T scale)
-      {
-        dolfinx::fem::set_bc<T>(std::span(b.data(), b.size()), bcs,
-                                std::span(x0.data(), x0.shape(0)), scale);
-      },
-      nb::arg("b").noconvert(), nb::arg("bcs"), nb::arg("x0").noconvert(),
-      nb::arg("scale"));
-  m.def(
-      "set_bc",
-      [](nb::ndarray<T, nb::ndim<1>, nb::c_contig> b,
-         const std::vector<
-             std::shared_ptr<const dolfinx::fem::DirichletBC<T, U>>>& bcs,
-         T scale)
-      { dolfinx::fem::set_bc<T>(std::span(b.data(), b.size()), bcs, scale); },
-      nb::arg("b").noconvert(), nb::arg("bcs"), nb::arg("scale"));
 }
 
 } // namespace
diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp
index 4725847162a..f28530b5e03 100644
--- a/python/dolfinx/wrappers/fem.cpp
+++ b/python/dolfinx/wrappers/fem.cpp
@@ -35,6 +35,7 @@
 #include <nanobind/stl/complex.h>
 #include <nanobind/stl/function.h>
 #include <nanobind/stl/map.h>
+#include <nanobind/stl/optional.h>
 #include <nanobind/stl/pair.h>
 #include <nanobind/stl/set.h>
 #include <nanobind/stl/shared_ptr.h>
@@ -109,7 +110,8 @@ void declare_function_space(nb::module_& m, std::string type)
             "__init__",
             [](dolfinx::fem::FiniteElement<T>* self,
                basix::FiniteElement<T>& element, std::size_t block_size,
-               bool symmetric) {
+               bool symmetric)
+            {
               new (self) dolfinx::fem::FiniteElement<T>(element, block_size,
                                                         symmetric);
             },
@@ -364,6 +366,24 @@ void declare_objects(nb::module_& m, const std::string& type)
                                   dofs.data(), {dofs.size()}, nb::handle()),
                               owned);
            })
+      .def(
+          "set",
+          [](const dolfinx::fem::DirichletBC<T, U>& self,
+             nb::ndarray<T, nb::ndim<1>, nb::c_contig> b,
+             std::optional<nb::ndarray<const T, nb::ndim<1>, nb::c_contig>> x0,
+             T alpha)
+          {
+            auto _b = std::span(b.data(), b.size());
+            if (x0.has_value())
+            {
+              self.set(_b, std::span(x0.value().data(), x0.value().shape(0)),
+                       alpha);
+            }
+            else
+              self.set(_b, std::nullopt, alpha);
+          },
+          nb::arg("b").noconvert(), nb::arg("x0").noconvert().none(),
+          nb::arg("alpha"))
       .def_prop_ro("function_space",
                    &dolfinx::fem::DirichletBC<T, U>::function_space)
       .def_prop_ro("value", &dolfinx::fem::DirichletBC<T, U>::value);
diff --git a/python/test/unit/fem/test_assemble_domains.py b/python/test/unit/fem/test_assemble_domains.py
index d454203ca9c..839fbcc31c4 100644
--- a/python/test/unit/fem/test_assemble_domains.py
+++ b/python/test/unit/fem/test_assemble_domains.py
@@ -92,13 +92,13 @@ def test_assembly_dx_domains(mode, meshtags_factory):
 
     fem.apply_lifting(b.array, [a], [[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b.array, [bc])
+    bc.set(b.array)
 
     L2 = form(ufl.inner(w, v) * dx)
     b2 = fem.assemble_vector(L2)
     fem.apply_lifting(b2.array, [a], [[bc]])
     b2.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b2.array, [bc])
+    bc.set(b2.array)
     assert np.allclose(b.array, b2.array)
 
     # Assemble scalar
@@ -181,13 +181,13 @@ def right(x):
 
     fem.apply_lifting(b.array, [a], [[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b.array, [bc])
+    bc.set(b.array)
 
     L2 = form(ufl.inner(w, v) * ds)
     b2 = fem.assemble_vector(L2)
     fem.apply_lifting(b2.array, [a2], [[bc]])
     b2.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b2.array, [bc])
+    bc.set(b2.array)
     assert np.allclose(b.array, b2.array)
 
     # Assemble scalar
diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py
index 95cafb12d0a..885661f60ee 100644
--- a/python/test/unit/fem/test_assemble_submesh.py
+++ b/python/test/unit/fem/test_assemble_submesh.py
@@ -58,7 +58,7 @@ def assemble(mesh, space, k):
     b = fem.assemble_vector(L)
     fem.apply_lifting(b.array, [a], bcs=[[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    fem.set_bc(b.array, [bc])
+    bc.set(b.array)
     s = mesh.comm.allreduce(
         fem.assemble_scalar(fem.form(ufl.inner(c * f, f) * (dx + ds))), op=MPI.SUM
     )
diff --git a/python/test/unit/fem/test_bcs.py b/python/test/unit/fem/test_bcs.py
index 2c57335a00d..fc658430733 100644
--- a/python/test/unit/fem/test_bcs.py
+++ b/python/test/unit/fem/test_bcs.py
@@ -25,7 +25,6 @@
     functionspace,
     locate_dofs_geometrical,
     locate_dofs_topological,
-    set_bc,
 )
 from dolfinx.mesh import CellType, create_unit_cube, create_unit_square, exterior_facet_indices
 from ufl import dx, inner
@@ -103,7 +102,8 @@ def test_overlapping_bcs():
     assemble_vector(b.array, L)
     apply_lifting(b.array, [a], [bcs])
     b.scatter_reverse(la.InsertMode.add)
-    set_bc(b.array, bcs)
+    for bc in bcs:
+        bc.set(b.array)
     b.scatter_forward()
 
     if len(dof_corner) > 0:
@@ -176,10 +176,10 @@ def test_constant_bc(mesh_factory):
     bc_c = dirichletbc(c, boundary_dofs, V)
 
     u_f = Function(V)
-    set_bc(u_f.x.array, [bc_f])
+    bc_f.set(u_f.x.array)
 
     u_c = Function(V)
-    set_bc(u_c.x.array, [bc_c])
+    bc_c.set(u_c.x.array)
     assert np.allclose(u_f.x.array, u_c.x.array)
 
 
@@ -217,14 +217,15 @@ def test_vector_constant_bc(mesh_factory):
         u_bcs[i].x.array[:] = c[i]
         bcs_f.append(dirichletbc(u_bcs[i], boundary_dofs[i], V.sub(i)))
     u_f = Function(V)
-    set_bc(u_f.x.array, bcs_f)
+    for bc in bcs_f:
+        bc.set(u_f.x.array)
 
     # Set using constant
     boundary_dofs = locate_dofs_topological(V, tdim - 1, boundary_facets)
     bc_c = dirichletbc(c, boundary_dofs, V)
     u_c = Function(V)
     u_c.x.array[:] = 0.0
-    set_bc(u_c.x.array, [bc_c])
+    bc_c.set(u_c.x.array)
 
     assert np.allclose(u_f.x.array, u_c.x.array)
 
@@ -262,9 +263,9 @@ def test_sub_constant_bc(mesh_factory):
         bc_c = dirichletbc(c, boundary_dofs, V.sub(i))
 
         u_f = Function(V)
-        set_bc(u_f.x.array, [bc_fi])
+        bc_fi.set(u_f.x.array)
         u_c = Function(V)
-        set_bc(u_c.x.array, [bc_c])
+        bc_c.set(u_c.x.array)
         assert np.allclose(u_f.x.array, u_c.x.array)
 
 
@@ -304,7 +305,7 @@ def test_mixed_constant_bc(mesh_factory):
         # Apply BC to scalar component of a mixed space using a Constant
         dofs = locate_dofs_topological(W.sub(i), tdim - 1, boundary_facets)
         bc = dirichletbc(c, dofs, W.sub(i))
-        set_bc(u.x.array, [bc])
+        bc.set(u.x.array)
 
         # Apply BC to scalar component of a mixed space using a Function
         ubc = u.sub(i).collapse()
@@ -313,7 +314,7 @@ def test_mixed_constant_bc(mesh_factory):
             (W.sub(i), ubc.function_space), tdim - 1, boundary_facets
         )
         bc_func = dirichletbc(ubc, dofs_both, W.sub(i))
-        set_bc(u_func.x.array, [bc_func])
+        bc_func.set(u_func.x.array)
 
         # Check that both approaches yield the same vector
         assert np.allclose(u.x.array, u_func.x.array)
@@ -344,7 +345,7 @@ def test_mixed_blocked_constant():
     c0 = default_scalar_type(3)
     dofs0 = locate_dofs_topological(W.sub(0), tdim - 1, boundary_facets)
     bc0 = dirichletbc(c0, dofs0, W.sub(0))
-    set_bc(u.x.array, [bc0])
+    bc0.set(u.x.array)
 
     # Apply BC to scalar component of a mixed space using a Function
     ubc = u.sub(0).collapse()
@@ -352,7 +353,7 @@ def test_mixed_blocked_constant():
     dofs_both = locate_dofs_topological((W.sub(0), ubc.function_space), tdim - 1, boundary_facets)
     bc_func = dirichletbc(ubc, dofs_both, W.sub(0))
     u_func = Function(W)
-    set_bc(u_func.x.array, [bc_func])
+    bc_func.set(u_func.x.array)
     assert np.allclose(u.x.array, u_func.x.array)
 
     # Check that vector space throws error
diff --git a/python/test/unit/fem/test_custom_basix_element.py b/python/test/unit/fem/test_custom_basix_element.py
index 8135db44522..a371b82ef35 100644
--- a/python/test/unit/fem/test_custom_basix_element.py
+++ b/python/test/unit/fem/test_custom_basix_element.py
@@ -23,7 +23,6 @@
     form,
     functionspace,
     locate_dofs_topological,
-    set_bc,
 )
 from dolfinx.mesh import CellType, create_unit_cube, create_unit_square, exterior_facet_indices
 from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, grad, inner
@@ -66,7 +65,7 @@ def run_scalar_test(V, degree, dtype, cg_solver, rtol=None):
     b = assemble_vector(L)
     apply_lifting(b.array, [a], bcs=[[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    set_bc(b.array, [bc])
+    bc.set(b.array)
 
     a = form(a, dtype=dtype)
     A = assemble_matrix(a, bcs=[bc])
diff --git a/python/test/unit/fem/test_fem_pipeline.py b/python/test/unit/fem/test_fem_pipeline.py
index 5179daefb30..d0501b2554e 100644
--- a/python/test/unit/fem/test_fem_pipeline.py
+++ b/python/test/unit/fem/test_fem_pipeline.py
@@ -26,7 +26,6 @@
     form,
     functionspace,
     locate_dofs_topological,
-    set_bc,
 )
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import (
@@ -93,7 +92,7 @@ def run_scalar_test(mesh, V, degree, cg_solver):
     b = assemble_vector(L)
     apply_lifting(b.array, [a], bcs=[[bc]])
     b.scatter_reverse(la.InsertMode.add)
-    set_bc(b.array, [bc])
+    bc.set(b.array)
 
     a = form(a, dtype=dtype)
     A = assemble_matrix(a, bcs=[bc])
@@ -386,7 +385,8 @@ def b(tau_S, v):
     b = assemble_vector(L)
     apply_lifting(b.array, [a], bcs=[bcs])
     b.scatter_reverse(la.InsertMode.add)
-    set_bc(b.array, bcs)
+    for bc in bcs:
+        bc.set(b.array)
 
     x_h = Function(V, dtype=dtype)
     x_h.x.array[:] = scipy.sparse.linalg.spsolve(A.to_scipy(), b.array)
diff --git a/python/test/unit/fem/test_petsc_nonlinear_assembler.py b/python/test/unit/fem/test_petsc_nonlinear_assembler.py
index 5f7bb7b3e92..5e79943feff 100644
--- a/python/test/unit/fem/test_petsc_nonlinear_assembler.py
+++ b/python/test/unit/fem/test_petsc_nonlinear_assembler.py
@@ -61,7 +61,7 @@ def F_mono(self, snes, x, F):
         with F.localForm() as f_local:
             f_local.set(0.0)
         assemble_vector(F, self.L)
-        apply_lifting(F, [self.a], bcs=[self.bcs], x0=[x], scale=-1.0)
+        apply_lifting(F, [self.a], bcs=[self.bcs], x0=[x], alpha=-1.0)
         F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
         set_bc(F, self.bcs, x, -1.0)
 
@@ -97,7 +97,7 @@ def F_block(self, snes, x, F):
             )
             offset += size_local
 
-        assemble_vector_block(F, self.L, self.a, bcs=self.bcs, x0=x, scale=-1.0)
+        assemble_vector_block(F, self.L, self.a, bcs=self.bcs, x0=x, alpha=-1.0)
 
     def J_block(self, snes, x, J, P):
         from dolfinx.fem.petsc import assemble_matrix_block
@@ -130,7 +130,7 @@ def F_nest(self, snes, x, F):
             with F_sub.localForm() as F_sub_local:
                 F_sub_local.set(0.0)
             assemble_vector(F_sub, L)
-            apply_lifting(F_sub, a, bcs=bcs1, x0=x, scale=-1.0)
+            apply_lifting(F_sub, a, bcs=bcs1, x0=x, alpha=-1.0)
             F_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
 
         # Set bc value in RHS
@@ -241,7 +241,7 @@ def blocked():
 
             # Ghosts are updated inside assemble_vector_block
             A = assemble_matrix_block(a_block, bcs=[bc])
-            b = assemble_vector_block(L_block, a_block, bcs=[bc], x0=x, scale=-1.0)
+            b = assemble_vector_block(L_block, a_block, bcs=[bc], x0=x, alpha=-1.0)
             A.assemble()
             assert A.getType() != "nest"
             Anorm = A.norm()
@@ -265,12 +265,12 @@ def nested():
 
             A = assemble_matrix_nest(a_block, bcs=[bc])
             b = assemble_vector_nest(L_block)
-            apply_lifting_nest(b, a_block, bcs=[bc], x0=x, scale=-1.0)
+            apply_lifting_nest(b, a_block, bcs=[bc], x0=x, alpha=-1.0)
             for b_sub in b.getNestSubVecs():
                 b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
             bcs0 = bcs_by_block([L.function_spaces[0] for L in L_block], [bc])
 
-            set_bc_nest(b, bcs0, x, scale=-1.0)
+            set_bc_nest(b, bcs0, x, alpha=-1.0)
             A.assemble()
             assert A.getType() == "nest"
             Anorm = nest_matrix_norm(A)
@@ -308,9 +308,9 @@ def monolithic():
             A = assemble_matrix(J, bcs=[bc])
             A.assemble()
             b = assemble_vector(F)
-            apply_lifting(b, [J], bcs=[[bc]], x0=[U.x.petsc_vec], scale=-1.0)
+            apply_lifting(b, [J], bcs=[[bc]], x0=[U.x.petsc_vec], alpha=-1.0)
             b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
-            set_bc(b, [bc], x0=U.x.petsc_vec, scale=-1.0)
+            set_bc(b, [bc], x0=U.x.petsc_vec, alpha=-1.0)
             assert A.getType() != "nest"
             Anorm = A.norm()
             bnorm = b.norm()
diff --git a/python/test/unit/nls/test_newton.py b/python/test/unit/nls/test_newton.py
index c458dbb6fd7..ba56dbcdf1a 100644
--- a/python/test/unit/nls/test_newton.py
+++ b/python/test/unit/nls/test_newton.py
@@ -42,7 +42,7 @@ def F(self, x, b):
         with b.localForm() as b_local:
             b_local.set(0.0)
         assemble_vector(b, self.L)
-        apply_lifting(b, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
+        apply_lifting(b, [self.a], bcs=[[self.bc]], x0=[x], alpha=-1.0)
         b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
         set_bc(b, [self.bc], x, -1.0)
 
@@ -88,7 +88,7 @@ def F(self, snes, x, F):
         with F.localForm() as f_local:
             f_local.set(0.0)
         assemble_vector(F, self.L)
-        apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
+        apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], alpha=-1.0)
         F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
         set_bc(F, [self.bc], x, -1.0)