Merge branch 'main' into main

Author: jorgensd

.github/workflows/sonarcloud.yml

@@ -17,7 +17,7 @@ jobs:
     container: fenicsproject/test-env:current-mpich
     env:
       SONAR_SCANNER_VERSION:
-        5.0.1.3006 # Find the latest version at:
+        6.1.0.4477 # Find the latest version at:
         # https://github.com/SonarSource/sonar-scanner-cli/tags
       SONAR_SERVER_URL: "https://sonarcloud.io"
       BUILD_WRAPPER_OUT_DIR: build_wrapper_output_directory # Directory where build-wrapper output will be placed
@@ -43,12 +43,12 @@ jobs:
           restore-keys: ${{ runner.os }}-sonar
       - name: Download and set up sonar-scanner
         env:
-          SONAR_SCANNER_DOWNLOAD_URL: https://binaries.sonarsource.com/Distribution/sonar-scanner-cli/sonar-scanner-cli-${{ env.SONAR_SCANNER_VERSION }}-linux.zip
+          SONAR_SCANNER_DOWNLOAD_URL: https://binaries.sonarsource.com/Distribution/sonar-scanner-cli/sonar-scanner-cli-${{ env.SONAR_SCANNER_VERSION }}-linux-x64.zip
         run: |
           mkdir -p $HOME/.sonar
           wget -O $HOME/.sonar/sonar-scanner.zip ${{ env.SONAR_SCANNER_DOWNLOAD_URL }}
           unzip -o $HOME/.sonar/sonar-scanner.zip -d $HOME/.sonar/
-          echo "$HOME/.sonar/sonar-scanner-${{ env.SONAR_SCANNER_VERSION }}-linux/bin" >> $GITHUB_PATH
+          echo "$HOME/.sonar/sonar-scanner-${{ env.SONAR_SCANNER_VERSION }}-linux-x64/bin" >> $GITHUB_PATH
       - name: Download and set up build-wrapper
         env:
           BUILD_WRAPPER_DOWNLOAD_URL: ${{ env.SONAR_SERVER_URL }}/static/cpp/build-wrapper-linux-x86.zip
@@ -71,4 +71,4 @@ jobs:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           SONAR_TOKEN: ${{ secrets.SONAR_TOKEN }}
         run: |
-          sonar-scanner --define sonar.host.url="${{ env.SONAR_SERVER_URL }}" --define sonar.cfamily.build-wrapper-output="${{ env.BUILD_WRAPPER_OUT_DIR }}"
+          sonar-scanner --define sonar.host.url="${{ env.SONAR_SERVER_URL }}" --define sonar.cfamily.compile-commands="${{ env.BUILD_WRAPPER_OUT_DIR }}/compile_commands.json"

cpp/demo/interpolation-io/main.cpp

@@ -29,8 +29,6 @@ using namespace dolfinx;
 /// and outputs the finite element function to a VTX file for
 /// visualisation.
 ///
-/// Also shows how to create a finite element using Basix.
-///
 /// @tparam T Scalar type of the finite element function.
 /// @tparam U Float type for the finite element basis and the mesh.
 /// @param mesh Mesh.
@@ -190,8 +188,8 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
 #endif
 }
 
-/// @brief This program shows how to create finite element spaces without FFCx
-/// generated code.
+/// @brief This program shows how to interpolate functions into different types
+/// of finite element spaces and output the result to file for visualisation.
 int main(int argc, char* argv[])
 {
   dolfinx::init_logging(argc, argv);

cpp/dolfinx/common/IndexMap.cpp

@@ -19,7 +19,6 @@ using namespace dolfinx::common;
 
 namespace
 {
-
 /// @brief Given source ranks (ranks that own indices ghosted by the
 /// calling rank), compute ranks that ghost indices owned by the calling
 /// rank.
@@ -156,6 +155,7 @@ communicate_ghosts_to_owners(MPI_Comm comm, std::span<const int> src,
 /// Given an index map and a subset of local indices (can be owned or
 /// ghost but must be unique and sorted), compute the owned, ghost and
 /// ghost owners in the submap.
+///
 /// @param[in] imap An index map.
 /// @param[in] indices List of entity indices (indices local to the
 /// process).
@@ -175,9 +175,6 @@ compute_submap_indices(const IndexMap& imap,
                        std::span<const std::int32_t> indices,
                        IndexMapOrder order, bool allow_owner_change)
 {
-  std::span<const int> src = imap.src();
-  std::span<const int> dest = imap.dest();
-
   // Create lookup array to determine if an index is in the sub-map
   std::vector<std::uint8_t> is_in_submap(imap.size_local() + imap.num_ghosts(),
                                          0);
@@ -210,7 +207,9 @@ compute_submap_indices(const IndexMap& imap,
     // processes that can own them in the submap.
     std::vector<std::pair<std::int64_t, int>> global_idx_to_possible_owner;
     const std::array local_range = imap.local_range();
+
     // Loop through the received indices
+    std::span<const int> dest = imap.dest();
     for (std::size_t i = 0; i < recv_disp.size() - 1; ++i)
     {
       for (int j = recv_disp[i]; j < recv_disp[i + 1]; ++j)
@@ -249,16 +248,16 @@ compute_submap_indices(const IndexMap& imap,
       // load balancing, though the impact is probably only very small
       auto it = std::ranges::lower_bound(global_idx_to_possible_owner, idx,
                                          std::ranges::less(),
-                                         [](auto& e) { return e.first; });
+                                         [](auto e) { return e.first; });
       assert(it != global_idx_to_possible_owner.end() and it->first == idx);
       send_owners.push_back(it->second);
     }
 
     // Create neighbourhood comm (owner -> ghost)
     MPI_Comm comm1;
     int ierr = MPI_Dist_graph_create_adjacent(
-        imap.comm(), src.size(), src.data(), MPI_UNWEIGHTED, dest.size(),
-        dest.data(), MPI_UNWEIGHTED, MPI_INFO_NULL, false, &comm1);
+        imap.comm(), imap.src().size(), imap.src().data(), MPI_UNWEIGHTED,
+        dest.size(), dest.data(), MPI_UNWEIGHTED, MPI_INFO_NULL, false, &comm1);
     dolfinx::MPI::check_error(imap.comm(), ierr);
 
     // Send the data
@@ -350,7 +349,7 @@ compute_submap_indices(const IndexMap& imap,
   }
 
   // Compute submap destination ranks
-  // FIXME Remove call to NBX
+  // FIXME: Remove call to NBX
   std::vector<int> submap_dest
       = dolfinx::MPI::compute_graph_edges_nbx(imap.comm(), submap_src);
   std::ranges::sort(submap_dest);
@@ -379,7 +378,7 @@ compute_submap_ghost_indices(std::span<const int> submap_src,
                              std::span<const std::int32_t> submap_owned,
                              std::span<const std::int64_t> submap_ghosts_global,
                              std::span<const std::int32_t> submap_ghost_owners,
-                             int submap_offset, const IndexMap& imap)
+                             std::int64_t submap_offset, const IndexMap& imap)
 {
   // --- Step 1 ---: Send global ghost indices (w.r.t. original imap) to
   // owning rank
@@ -397,15 +396,15 @@ compute_submap_ghost_indices(std::span<const int> submap_src,
   std::vector<std::int64_t> send_gidx;
   {
     send_gidx.reserve(recv_indices.size());
-    // NOTE: Received indices are owned by this process in the submap, but not
-    // necessarily in the original imap, so we must use global_to_local to
-    // convert rather than subtracting local_range[0]
-    // TODO Convert recv_indices or submap_owned?
-    std::vector<int32_t> recv_indices_local(recv_indices.size());
+    // NOTE: Received indices are owned by this process in the submap,
+    // but not necessarily in the original imap, so we must use
+    // global_to_local to convert rather than subtracting local_range[0]
+    // TODO: Convert recv_indices or submap_owned?
+    std::vector<std::int32_t> recv_indices_local(recv_indices.size());
     imap.global_to_local(recv_indices, recv_indices_local);
 
     // Compute submap global index
-    for (auto idx : recv_indices_local)
+    for (std::int32_t idx : recv_indices_local)
     {
       // Could avoid search by creating look-up array
       auto it = std::ranges::lower_bound(submap_owned, idx);
@@ -415,7 +414,8 @@ compute_submap_ghost_indices(std::span<const int> submap_src,
     }
   }
 
-  // --- Step 3 ---: Send submap global indices to process that ghost them
+  // --- Step 3 ---: Send submap global indices to process that ghost
+  // them
 
   std::vector<std::int64_t> recv_gidx(send_disp.back());
   {

cpp/dolfinx/fem/Expression.h

@@ -199,7 +199,7 @@ class Expression
       num_argument_dofs
           = _argument_function_space->dofmap()->element_dof_layout().num_dofs();
       auto element = _argument_function_space->element();
-
+      num_argument_dofs *= _argument_function_space->dofmap()->bs();
       assert(element);
       if (element->needs_dof_transformations())
       {
@@ -244,7 +244,6 @@ class Expression
       std::ranges::fill(values_local, 0);
       _fn(values_local.data(), coeff_cell, constant_data.data(),
           coord_dofs.data(), entity_index, nullptr);
-
       post_dof_transform(values_local, cell_info, e, size0);
       for (std::size_t j = 0; j < values_local.size(); ++j)
         values[e * vshape[1] + j] = values_local[j];

cpp/dolfinx/fem/interpolate.h

@@ -134,12 +134,12 @@ using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
 
 /// @brief Scatter data into non-contiguous memory.
 ///
-/// Scatter blocked data `send_values` to its corresponding `src_rank` and
-/// insert the data into `recv_values`. The insert location in
+/// Scatter blocked data `send_values` to its corresponding `src_rank`
+/// and insert the data into `recv_values`. The insert location in
 /// `recv_values` is determined by `dest_ranks`. If the j-th dest rank
 /// is -1, then `recv_values[j*block_size:(j+1)*block_size]) = 0`.
 ///
-/// @param[in] comm The MPI communicator
+/// @param[in] comm The MPI communicator.
 /// @param[in] src_ranks Rank owning the values of each row in
 /// `send_values`.
 /// @param[in] dest_ranks List of ranks receiving data. Size of array is
@@ -194,9 +194,11 @@ void scatter_values(MPI_Comm comm, std::span<const std::int32_t> src_ranks,
   std::vector<std::int32_t> recv_offsets(in_ranks.size() + 1, 0);
   {
     // Build map from parent to neighborhood communicator ranks
-    std::map<std::int32_t, std::int32_t> rank_to_neighbor;
+    std::vector<std::pair<std::int32_t, std::int32_t>> rank_to_neighbor;
+    rank_to_neighbor.reserve(in_ranks.size());
     for (std::size_t i = 0; i < in_ranks.size(); i++)
-      rank_to_neighbor[in_ranks[i]] = i;
+      rank_to_neighbor.push_back({in_ranks[i], i});
+    std::ranges::sort(rank_to_neighbor);
 
     // Compute receive sizes
     std::ranges::for_each(
@@ -205,8 +207,11 @@ void scatter_values(MPI_Comm comm, std::span<const std::int32_t> src_ranks,
         {
           if (rank >= 0)
           {
-            const int neighbor = rank_to_neighbor[rank];
-            recv_sizes[neighbor] += block_size;
+            auto it = std::ranges::lower_bound(rank_to_neighbor, rank,
+                                               std::ranges::less(),
+                                               [](auto e) { return e.first; });
+            assert(it != rank_to_neighbor.end() and it->first == rank);
+            recv_sizes[it->second] += block_size;
           }
         });
 
@@ -221,26 +226,42 @@ void scatter_values(MPI_Comm comm, std::span<const std::int32_t> src_ranks,
     {
       if (const std::int32_t rank = dest_ranks[i]; rank >= 0)
       {
-        const int neighbor = rank_to_neighbor[rank];
-        int insert_pos = recv_offsets[neighbor] + recv_counter[neighbor];
+        auto it = std::ranges::lower_bound(rank_to_neighbor, rank,
+                                           std::ranges::less(),
+                                           [](auto e) { return e.first; });
+        assert(it != rank_to_neighbor.end() and it->first == rank);
+        int insert_pos = recv_offsets[it->second] + recv_counter[it->second];
         comm_to_output[insert_pos / block_size] = i * block_size;
-        recv_counter[neighbor] += block_size;
+        recv_counter[it->second] += block_size;
       }
     }
   }
 
   std::vector<std::int32_t> send_sizes(out_ranks.size());
   send_sizes.reserve(1);
   {
-    // Compute map from parent mpi rank to neigbor rank for outgoing data
-    std::map<std::int32_t, std::int32_t> rank_to_neighbor;
+    // Compute map from parent MPI rank to neighbor rank for outgoing
+    // data. `out_ranks` is sorted, so rank_to_neighbor will be sorted
+    // too.
+    std::vector<std::pair<std::int32_t, std::int32_t>> rank_to_neighbor;
+    rank_to_neighbor.reserve(out_ranks.size());
     for (std::size_t i = 0; i < out_ranks.size(); i++)
-      rank_to_neighbor[out_ranks[i]] = i;
+      rank_to_neighbor.push_back({out_ranks[i], i});
 
-    // Compute send sizes
+    // Compute send sizes. As `src_ranks` is sorted, we can move 'start'
+    // in search forward.
+    auto start = rank_to_neighbor.begin();
     std::ranges::for_each(
-        src_ranks, [&rank_to_neighbor, &send_sizes, block_size](auto rank)
-        { send_sizes[rank_to_neighbor[rank]] += block_size; });
+        src_ranks,
+        [&rank_to_neighbor, &send_sizes, block_size, &start](auto rank)
+        {
+          auto it = std::ranges::lower_bound(start, rank_to_neighbor.end(),
+                                             rank, std::ranges::less(),
+                                             [](auto e) { return e.first; });
+          assert(it != rank_to_neighbor.end() and it->first == rank);
+          send_sizes[it->second] += block_size;
+          start = it;
+        });
   }
 
   // Compute sending offsets
@@ -257,7 +278,8 @@ void scatter_values(MPI_Comm comm, std::span<const std::int32_t> src_ranks,
                          dolfinx::MPI::mpi_type<T>(), reverse_comm);
   MPI_Comm_free(&reverse_comm);
 
-  // Insert values received from neighborhood communicator in output span
+  // Insert values received from neighborhood communicator in output
+  // span
   std::ranges::fill(recv_values, T(0));
   for (std::size_t i = 0; i < comm_to_output.size(); i++)
   {
@@ -1117,10 +1139,11 @@ void interpolate(Function<T, U>& u, const Function<T, U>& v,
   assert(element_u);
   const std::size_t value_size = u.function_space()->value_size();
 
-  auto& dest_ranks = interpolation_data.src_owner;
-  auto& src_ranks = interpolation_data.dest_owners;
-  auto& recv_points = interpolation_data.dest_points;
-  auto& evaluation_cells = interpolation_data.dest_cells;
+  const std::vector<int>& dest_ranks = interpolation_data.src_owner;
+  const std::vector<int>& src_ranks = interpolation_data.dest_owners;
+  const std::vector<U>& recv_points = interpolation_data.dest_points;
+  const std::vector<std::int32_t>& evaluation_cells
+      = interpolation_data.dest_cells;
 
   // Evaluate the interpolating function where possible
   std::vector<T> send_values(recv_points.size() / 3 * value_size);

cpp/dolfinx/mesh/Geometry.h

@@ -113,7 +113,7 @@ class Geometry
   /// Move Assignment
   Geometry& operator=(Geometry&&) = default;
 
-  /// Return Euclidean dimension of coordinate system
+  /// Return dimension of the Euclidean coordinate system
   int dim() const { return _dim; }
 
   /// @brief  DofMap for the geometry

cpp/dolfinx/mesh/Topology.h

@@ -151,16 +151,18 @@ class Topology
   /// @brief Returns the permutation information
   const std::vector<std::uint32_t>& get_cell_permutation_info() const;
 
-  /// @brief Get the permutation number to apply to a facet.
+  /// @brief Get the numbers that encode the number of permutations to apply to
+  /// facets.
   ///
-  /// The permutations are numbered so that:
+  /// The permutations are encoded so that:
   ///
   ///   - `n % 2` gives the number of reflections to apply
   ///   - `n // 2` gives the number of rotations to apply
   ///
-  /// Each column of the returned array represents a cell, and each row
-  /// a facet of that cell.
-  /// @return The permutation number
+  /// The data is stored in a flattened 2D array, so that `data[cell_index *
+  /// facets_per_cell + facet_index]` contains the facet with index
+  /// `facet_index` of the cell with index `cell_index`.
+  /// @return The encoded permutation info
   /// @note An exception is raised if the permutations have not been
   /// computed
   const std::vector<std::uint8_t>& get_facet_permutations() const;