Demo updates for LU solver bugs (#3499)

* Docker update

* Type fixes in pyvista demo

* Change solver type

* More solver switches

* Demo update

* Run all tests

* Test update

* Try mumps option

* Demo update

* Demo update

* Update

* Update

* mypy fix

* Update

* Try superlu

* Lint

* Demo update

* Demo update

* Updates

* Update

* Update

* Update

* Demo update

* Update

* Update

* Lint

* mypy fixex

* Tidy up

* Lint

* Formatting fixes

* Logic fix

* Formatting fixes

.github/workflows/ccpp.yml

@@ -33,7 +33,7 @@ jobs:
     name: Lint
     steps:
       - uses: actions/checkout@v4
-      
+
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
@@ -85,7 +85,7 @@ jobs:
     name: Build and test
     steps:
       - uses: actions/checkout@v4
-      
+
       - name: Install dependencies
         run: |
           sudo apt-get update
@@ -209,8 +209,8 @@ jobs:
           ctest -V -R demo -R serial
           ctest -V -R demo -R mpi_2
 
-      - name: Install Python build dependencies 
-        run: pip install -r python/build-requirements.txt 
+      - name: Install Python build dependencies
+        run: pip install -r python/build-requirements.txt
 
       - name: Build Python interface
         run: |
@@ -282,15 +282,15 @@ jobs:
         uses: actions/upload-artifact@v4
         with:
           name: docs-cpp-doxygen
-          path: cpp/doc/html 
-          retention-days: 2 
+          path: cpp/doc/html
+          retention-days: 2
       - name: Upload C++ Sphinx documentation artifact
         uses: actions/upload-artifact@v4
         with:
           name: docs-cpp-sphinx
-          path: cpp/doc/build/html 
-          retention-days: 2 
-      
+          path: cpp/doc/build/html
+          retention-days: 2
+
       - name: Build Python interface documentation
         run: |
           cd python/doc
@@ -299,8 +299,8 @@ jobs:
         uses: actions/upload-artifact@v4
         with:
           name: docs-python
-          path: python/doc/build/html 
-          retention-days: 2 
+          path: python/doc/build/html
+          retention-days: 2
 
       - name: Checkout FEniCS/docs
         if: ${{ github.repository == 'FEniCS/dolfinx' && ( github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/v') ) }}

cpp/dolfinx/common/Scatterer.h

@@ -145,7 +145,8 @@ class Scatterer
 
     // Scale sizes and displacements by block size
     {
-      auto rescale = [](auto& x, int bs) {
+      auto rescale = [](auto& x, int bs)
+      {
         std::ranges::transform(x, x.begin(), [bs](auto e) { return e *= bs; });
       };
       rescale(_sizes_local, bs);
@@ -459,8 +460,8 @@ class Scatterer
   /// @note This function is intended for advanced usage, and in
   /// particular when using CUDA/device-aware MPI.
   ///
-  /// @tparam T The data type to send
-  /// @tparam F The pack function
+  /// @tparam T Data type to send.
+  /// @tparam F Pack function.
   /// @param[in] remote_data Received data associated with the ghost
   /// indices. The order follows the order of the ghost indices in the
   /// IndexMap used to create the scatterer. The size equal to the
@@ -493,14 +494,15 @@ class Scatterer
                       type);
   }
 
-  /// @brief End the reverse scatter communication, and unpack the received
-  /// local buffer into local data.
+  /// @brief End the reverse scatter communication, and unpack the
+  /// received local buffer into local data.
   ///
   /// This function must be called after Scatterer::scatter_rev_begin.
   /// The buffers passed to Scatterer::scatter_rev_begin must not be
   /// modified until after the function has been called.
-  /// @param[in] local_buffer Working buffer. Same buffer should be used in
-  /// Scatterer::scatter_rev_begin.
+  ///
+  /// @param[in] local_buffer Working buffer. Same buffer should be used
+  /// in Scatterer::scatter_rev_begin.
   /// @param[out] local_data All data associated with owned indices.
   /// Size is `size_local()` from the IndexMap used to create the
   /// scatterer, multiplied by the block size. The data for each index
@@ -557,7 +559,7 @@ class Scatterer
   }
 
   /// @brief Size of buffer for local data (owned and shared) used in
-  /// forward and reverse communication
+  /// forward and reverse communication.
   /// @return The required buffer size
   std::int32_t local_buffer_size() const noexcept { return _local_inds.size(); }
 
@@ -569,9 +571,9 @@ class Scatterer
     return _remote_inds.size();
   }
 
-  /// Return a vector of local indices (owned) used to pack/unpack local data.
-  /// These indices are grouped by neighbor process (process for which an index
-  /// is a ghost).
+  /// Return a vector of local indices (owned) used to pack/unpack local
+  /// data. These indices are grouped by neighbor process (process for
+  /// which an index is a ghost).
   const std::vector<std::int32_t>& local_indices() const noexcept
   {
     return _local_inds;

cpp/dolfinx/io/xdmf_utils.cpp

@@ -403,7 +403,8 @@ xdmf_utils::distribute_entity_data(
     std::vector<std::pair<int, std::int64_t>> dest_to_index;
     std::ranges::transform(
         indices, std::back_inserter(dest_to_index),
-        [size, num_nodes](auto n) {
+        [size, num_nodes](auto n)
+        {
           return std::pair(dolfinx::MPI::index_owner(size, n, num_nodes), n);
         });
     std::ranges::sort(dest_to_index);

cpp/dolfinx/io/xdmf_utils.h

@@ -26,7 +26,6 @@
 
 namespace dolfinx
 {
-
 namespace fem
 {
 template <dolfinx::scalar T, std::floating_point U>

cpp/dolfinx/la/MatrixCSR.h

@@ -70,17 +70,18 @@ class MatrixCSR
   /// @brief Insertion functor for setting values in a matrix. It is
   /// typically used in finite element assembly functions.
   ///
-  /// Create a function to set values in a MatrixCSR. The function signature is
-  /// `int mat_set_fn(std::span<const std::int32_t rows, std::span<const
-  /// std::int32_t cols, std::span<const value_type> data)`. The rows
-  /// and columns use process local indexing, and the given rows and
-  /// columns must pre-exist in the sparsity pattern of the matrix.
-  /// Insertion into "ghost" rows (in the ghost region of the row
-  /// `IndexMap`) is permitted, so long as there are correct entries
-  /// in the sparsity pattern.
+  /// Create a function to set values in a MatrixCSR. The function
+  /// signature is `int mat_set_fn(std::span<const std::int32_t rows,
+  /// std::span<const std::int32_t cols, std::span<const value_type>
+  /// data)`. The rows and columns use process local indexing, and the
+  /// given rows and columns must pre-exist in the sparsity pattern of
+  /// the matrix. Insertion into "ghost" rows (in the ghost region of
+  /// the row `IndexMap`) is permitted, so long as there are correct
+  /// entries in the sparsity pattern.
   ///
-  /// @note Using rows or columns which are not in the sparsity will result
-  /// in undefined behaviour (or an assert failure in Debug mode).
+  /// @note Using rows or columns which are not in the sparsity will
+  /// result in undefined behaviour (or an assert failure in Debug
+  /// mode).
   ///
   /// @note Matrix block size may be (1, 1) or (BS0, BS1)
   /// @note Data block size may be (1, 1) or (BS0, BS1)
@@ -111,17 +112,18 @@ class MatrixCSR
   /// @brief Insertion functor for adding values to a matrix. It is
   /// typically used in finite element assembly functions.
   ///
-  /// Create a function to add values to a MatrixCSR. The function signature is
-  /// `int mat_add_fn(std::span<const std::int32_t rows, std::span<const
-  /// std::int32_t cols, std::span<const value_type> data)`. The rows
-  /// and columns use process local indexing, and the given rows and
-  /// columns must pre-exist in the sparsity pattern of the matrix.
-  /// Insertion into "ghost" rows (in the ghost region of the row
-  /// `IndexMap`) is permitted, so long as there are correct entries
-  /// in the sparsity pattern.
+  /// Create a function to add values to a MatrixCSR. The function
+  /// signature is `int mat_add_fn(std::span<const std::int32_t rows,
+  /// std::span<const std::int32_t cols, std::span<const value_type>
+  /// data)`. The rows and columns use process local indexing, and the
+  /// given rows and columns must pre-exist in the sparsity pattern of
+  /// the matrix. Insertion into "ghost" rows (in the ghost region of
+  /// the row `IndexMap`) is permitted, so long as there are correct
+  /// entries in the sparsity pattern.
   ///
-  /// @note Using rows or columns which are not in the sparsity will result
-  /// in undefined behaviour (or an assert failure in Debug mode).
+  /// @note Using rows or columns which are not in the sparsity will
+  /// result in undefined behaviour (or an assert failure in Debug
+  /// mode).
   ///
   /// @note Matrix block size may be (1, 1) or (BS0, BS1)
   /// @note Data block size may be (1, 1) or (BS0, BS1)
@@ -152,16 +154,14 @@ class MatrixCSR
   /// @brief Create a distributed matrix.
   ///
   /// The structure of the matrix depends entirely on the input
-  /// `SparsityPattern`, which must be finalized.
-  /// The matrix storage is distributed Compressed Sparse Row:
-  /// the matrix is distributed by row across processes, and on each
-  /// process, there is a list of column indices and matrix entries
-  /// for each row stored. This exactly matches the layout of the
-  /// `SparsityPattern`.
-  /// There is some overlap of matrix rows between processes to allow for
-  /// independent Finite Element assembly, after which, the ghost rows
-  /// should be sent to the row owning processes by calling
-  /// `scatter_rev()`.
+  /// `SparsityPattern`, which must be finalized. The matrix storage is
+  /// distributed Compressed Sparse Row: the matrix is distributed by
+  /// row across processes, and on each process, there is a list of
+  /// column indices and matrix entries for each row stored. This
+  /// exactly matches the layout of the `SparsityPattern`. There is some
+  /// overlap of matrix rows between processes to allow for independent
+  /// Finite Element assembly, after which, the ghost rows should be
+  /// sent to the row owning processes by calling `scatter_rev()`.
   ///
   /// @note The block size of the matrix is given by the block size of
   /// the input `SparsityPattern`.
@@ -192,8 +192,9 @@ class MatrixCSR
   /// @note All indices are local to the calling MPI rank and entries
   /// cannot be set in ghost rows.
   /// @note This should be called after `scatter_rev`. Using before
-  /// `scatter_rev` will set the values correctly, but incoming values may
-  /// get added to them during a subsequent reverse scatter operation.
+  /// `scatter_rev` will set the values correctly, but incoming values
+  /// may get added to them during a subsequent reverse scatter
+  /// operation.
   /// @tparam BS0 Data row block size
   /// @tparam BS1 Data column block size
   /// @param[in] x The `m` by `n` dense block of values (row-major) to
@@ -216,8 +217,8 @@ class MatrixCSR
     }
     else if (_bs[0] == 1 and _bs[1] == 1)
     {
-      // Set blocked data in a regular CSR matrix (_bs[0]=1, _bs[1]=1) with
-      // correct sparsity
+      // Set blocked data in a regular CSR matrix (_bs[0]=1, _bs[1]=1)
+      // with correct sparsity
       impl::insert_blocked_csr<BS0, BS1>(_data, _cols, _row_ptr, x, rows, cols,
                                          set_fn, num_rows);
     }
@@ -278,20 +279,23 @@ class MatrixCSR
   /// Number of local rows including ghost rows
   std::int32_t num_all_rows() const { return _row_ptr.size() - 1; }
 
-  /// Copy to a dense matrix
+  /// @brief Copy to a dense matrix.
   /// @note This function is typically used for debugging and not used
-  /// in production
+  /// in production.
   /// @note Ghost rows are also returned, and these can be truncated
   /// manually by using num_owned_rows() if required.
-  /// @note If the block size is greater than 1, the entries are expanded.
+  /// @note If the block size is greater than 1, the entries are
+  /// expanded.
   /// @return Dense copy of the part of the matrix on the calling rank.
   /// Storage is row-major.
   std::vector<value_type> to_dense() const;
 
   /// @brief Transfer ghost row data to the owning ranks accumulating
   /// received values on the owned rows, and zeroing any existing data
-  /// in ghost rows. This process is analogous to `scatter_rev` for `Vector`
-  /// except that the values are always accumulated on the owning process.
+  /// in ghost rows.
+  ///
+  /// This process is analogous to `scatter_rev` for `Vector` except
+  /// that the values are always accumulated on the owning process.
   void scatter_rev()
   {
     scatter_rev_begin();
@@ -303,8 +307,8 @@ class MatrixCSR
   /// @note Calls to this function must be followed by
   /// MatrixCSR::scatter_rev_end(). Between the two calls matrix values
   /// must not be changed.
-  /// @note This function does not change the matrix data. Data update only
-  /// occurs with `scatter_rev_end()`.
+  /// @note This function does not change the matrix data. Data update
+  /// only occurs with `scatter_rev_end()`.
   void scatter_rev_begin();
 
   /// @brief End transfer of ghost row data to owning ranks.
@@ -637,7 +641,7 @@ MatrixCSR<U, V, W, X>::to_dense() const
       for (int i0 = 0; i0 < _bs[0]; ++i0)
         for (int i1 = 0; i1 < _bs[1]; ++i1)
         {
-          std::array<std::int32_t, 1> local_col {_cols[j]};
+          std::array<std::int32_t, 1> local_col{_cols[j]};
           std::array<std::int64_t, 1> global_col{0};
           _index_maps[1]->local_to_global(local_col, global_col);
           A[(r * _bs[1] + i0) * ncols * _bs[0] + global_col[0] * _bs[1] + i1]

cpp/dolfinx/la/SparsityPattern.h

@@ -74,11 +74,11 @@ class SparsityPattern
   /// @brief Insert non-zero locations using local (process-wise)
   /// indices.
   ///
-  /// This routine inserts non-zero locations at the outer product of rows and
-  /// cols into the sparsity pattern, i.e. adds the matrix entries at
-  ///   A[row[i], col[j]] for all i, j.
+  /// This routine inserts non-zero locations at the outer product of
+  /// rows and cols into the sparsity pattern, i.e. adds the matrix
+  /// entries at `A[row[i], col[j]] for all i, j`.
   ///
-  /// @param[in] rows list of the local row indices 
+  /// @param[in] rows list of the local row indices
   /// @param[in] cols list of the local column indices
   void insert(std::span<const std::int32_t> rows,
               std::span<const std::int32_t> cols);
@@ -88,7 +88,8 @@ class SparsityPattern
   /// must exist in the row IndexMap.
   void insert_diagonal(std::span<const std::int32_t> rows);
 
-  /// @brief Finalize sparsity pattern and communicate off-process entries
+  /// @brief Finalize sparsity pattern and communicate off-process
+  /// entries
   void finalize();
 
   /// @brief Index map for given dimension dimension. Returns the index

cpp/dolfinx/mesh/Geometry.h

@@ -165,14 +165,14 @@ class Geometry
   /// @brief Access geometry degrees-of-freedom data (const version).
   ///
   /// @return The flattened row-major geometry data, where the shape is
-  /// (num_points, 3)
+  /// `(num_points, 3)`.
   std::span<const value_type> x() const { return _x; }
 
   /// @brief Access geometry degrees-of-freedom data (non-const
   /// version).
   ///
   /// @return The flattened row-major geometry data, where the shape is
-  /// (num_points, 3)
+  /// `(num_points, 3)`.
   std::span<value_type> x() { return _x; }
 
   /// @brief The element that describes the geometry map.
@@ -231,8 +231,8 @@ template <typename U, typename V, typename W>
 Geometry(std::shared_ptr<const common::IndexMap>, U,
          const std::vector<fem::CoordinateElement<
              typename std::remove_reference_t<typename V::value_type>>>&,
-         V, int,
-         W) -> Geometry<typename std::remove_cvref_t<typename V::value_type>>;
+         V, int, W)
+    -> Geometry<typename std::remove_cvref_t<typename V::value_type>>;
 /// @endcond
 
 /// @brief Build Geometry from input data.

cpp/dolfinx/mesh/Mesh.h

@@ -103,8 +103,8 @@ class Mesh
 /// @cond
 /// Template type deduction
 template <typename V>
-Mesh(MPI_Comm, std::shared_ptr<Topology>,
-     V) -> Mesh<typename std::remove_cvref_t<typename V::value_type>>;
+Mesh(MPI_Comm, std::shared_ptr<Topology>, V)
+    -> Mesh<typename std::remove_cvref_t<typename V::value_type>>;
 /// @endcond
 
 } // namespace dolfinx::mesh

cpp/dolfinx/mesh/topologycomputation.cpp

@@ -538,10 +538,9 @@ compute_entities_by_key_matching(
 
         std::vector<std::size_t> perm(global_vertices.size());
         std::iota(perm.begin(), perm.end(), 0);
-        std::ranges::sort(perm,
-                          [&global_vertices](std::size_t i0, std::size_t i1) {
-                            return global_vertices[i0] < global_vertices[i1];
-                          });
+        std::ranges::sort(
+            perm, [&global_vertices](std::size_t i0, std::size_t i1)
+            { return global_vertices[i0] < global_vertices[i1]; });
         // For quadrilaterals, the vertex opposite the lowest vertex should
         // be last
         if (entity_type == mesh::CellType::quadrilateral)