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

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

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

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

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

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

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

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