Merge branch 'jpdean/mixed_mesh_hackathon_2' of github.com:FEniCS/dol…

…finx into jpdean/mixed_mesh_hackathon_2

.github/workflows/ccpp.yml

@@ -78,7 +78,7 @@ jobs:
           find . -type f \( -name "*.cmake" -o -name "*.cmake.in" -o -name "CMakeLists.txt" \) | xargs cmake-format --check
 
   build:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-24.04
     needs: [fenicsx-refs, lint]
     env:
       OMPI_ALLOW_RUN_AS_ROOT: 1
@@ -93,7 +93,9 @@ jobs:
       - name: Install dependencies
         run: |
           sudo apt-get update
-          sudo apt-get install catch2 cmake g++ libboost-dev libhdf5-mpi-dev libparmetis-dev libpugixml-dev libspdlog-dev mpi-default-dev ninja-build pkg-config
+          sudo apt-get install catch2 cmake g++ libblas-dev libboost-dev libhdf5-mpi-dev \
+             liblapack-dev libparmetis-dev libpugixml-dev libspdlog-dev mpi-default-dev \
+             ninja-build pkg-config
       - name: Set up Python
         uses: actions/setup-python@v5
         with:

.github/workflows/macos.yml

@@ -20,7 +20,7 @@ jobs:
 
   mac-os-build:
     name: macOS Homebrew install and test
-    runs-on: macos-14
+    runs-on: macos-15
     needs: fenicsx-refs
     env:
       PETSC_ARCH: arch-darwin-c-opt

cpp/dolfinx/fem/CoordinateElement.h

@@ -78,7 +78,9 @@ class CoordinateElement
   /// compute. Use 0 for the basis functions only
   /// @param[in] num_points Number of points at which to evaluate the
   /// basis functions.
-  /// @return Shape of the array to be filled by `tabulate`.
+  /// @return Shape of the array to be filled by `tabulate`, where (0)
+  /// is derivative index, (1) is the point index, (2) is the basis
+  /// function index and (3) is the basis function component.
   std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
                                             std::size_t num_points) const;
 

cpp/dolfinx/fem/Form.h

@@ -277,9 +277,11 @@ class Form
     auto end = std::ranges::upper_bound(integrals, i, std::less<>{}, get_id);
 
     // Check that the kernel is valid and return it if so
-    if (start == integrals.end() || start->id != i
+    if (start == integrals.end() or start->id != i
         or std::distance(start, end) <= kernel_idx)
+    {
       throw std::runtime_error("No kernel for requested domain index.");
+    }
 
     return std::next(start, kernel_idx)->kernel;
   }
@@ -395,15 +397,16 @@ class Form
                                    int kernel_idx) const
   {
     const auto& integrals = _integrals[static_cast<std::size_t>(type)];
-
     auto get_id = [](const auto& a) { return a.id; };
     auto start = std::ranges::lower_bound(integrals, i, std::less<>{}, get_id);
     auto end = std::ranges::upper_bound(integrals, i, std::less<>{}, get_id);
 
     // Check that the kernel is valid and return it if so
-    if (start == integrals.end() || start->id != i
+    if (start == integrals.end() or start->id != i
         or std::distance(start, end) <= kernel_idx)
+    {
       throw std::runtime_error("No kernel for requested domain index.");
+    }
 
     return std::next(start, kernel_idx)->entities;
   }
@@ -414,8 +417,8 @@ class Form
   ///
   /// @param type Integral type.
   /// @param i Integral ID, i.e. the (sub)domain index.
-  /// @param kernel_idx Index of the kernel (we may have multiple kernels
-  /// for a given ID in mixed-topology meshes).
+  /// @param kernel_idx Index of the kernel (we may have multiple
+  /// kernels for a given ID in mixed-topology meshes).
   /// @param mesh The mesh the entities are numbered with respect to.
   /// @return List of active entities in `mesh` for the given integral.
   std::vector<std::int32_t> domain(IntegralType type, int i, int kernel_idx,
@@ -468,6 +471,7 @@ class Form
             // Get the facet index
             const std::int32_t facet
                 = c_to_f->links(entities[i])[entities[i + 1]];
+
             // Add cell and the local facet index
             mapped_entities.insert(mapped_entities.end(),
                                    {entity_map[facet], entities[i + 1]});
@@ -495,16 +499,17 @@ class Form
         }
         else if (codim == 1)
         {
-          // In this case, the entity maps take facets in (`_mesh`) to cells in
-          // `mesh`, so we need to get the facet number from the (cell,
-          // local_facet pair) first.
+          // In this case, the entity maps take facets in (`_mesh`) to
+          // cells in `mesh`, so we need to get the facet number from
+          // the (cell, local_facet pair) first.
           auto c_to_f = _mesh->topology()->connectivity(tdim, tdim - 1);
           assert(c_to_f);
           for (std::size_t i = 0; i < entities.size(); i += 2)
           {
             // Get the facet index
             const std::int32_t facet
                 = c_to_f->links(entities[i])[entities[i + 1]];
+
             // Add cell and the local facet index
             mapped_entities.insert(mapped_entities.end(),
                                    {entity_map[facet], entities[i + 1]});
@@ -597,5 +602,5 @@ class Form
   std::map<std::shared_ptr<const mesh::Mesh<geometry_type>>,
            std::vector<std::int32_t>>
       _entity_maps;
-}; // namespace dolfinx::fem
+};
 } // namespace dolfinx::fem

cpp/dolfinx/fem/FunctionSpace.h

@@ -344,8 +344,11 @@ class FunctionSpace
   std::shared_ptr<const FiniteElement<geometry_type>> element() const
   {
     if (_elements.size() > 1)
+    {
       throw std::runtime_error(
           "FunctionSpace has multiple elements, call `elements` instead.");
+    }
+
     return elements(0);
   }
 
@@ -360,8 +363,11 @@ class FunctionSpace
   std::shared_ptr<const DofMap> dofmap() const
   {
     if (_dofmaps.size() > 1)
+    {
       throw std::runtime_error(
           "FunctionSpace has multiple dofmaps, call `dofmaps` instead.");
+    }
+
     return dofmaps(0);
   }
 

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -510,11 +510,12 @@ void assemble_matrix(
   auto mesh1 = a.function_spaces().at(1)->mesh();
   assert(mesh1);
 
-  // TODO Mixed topology with exterior and interior facet integrals.
+  // TODO: Mixed topology with exterior and interior facet integrals.
+  //
   // NOTE: Can't just loop over cell types for interior facet integrals
   // because we have a kernel per combination of comparable cell types,
-  // rather than one per cell type. Also, we need the dofmaps for
-  // two different cell types at the same time.
+  // rather than one per cell type. Also, we need the dofmaps for two
+  // different cell types at the same time.
   const int num_cell_types = mesh->topology()->cell_types().size();
   for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
   {
@@ -583,8 +584,11 @@ void assemble_matrix(
     for (int i : a.integral_ids(IntegralType::exterior_facet))
     {
       if (num_cell_types > 1)
+      {
         throw std::runtime_error("Exterior facet integrals with mixed "
                                  "topology aren't supported yet");
+      }
+
       auto fn = a.kernel(IntegralType::exterior_facet, i);
       assert(fn);
       auto& [coeffs, cstride]
@@ -601,8 +605,11 @@ void assemble_matrix(
     for (int i : a.integral_ids(IntegralType::interior_facet))
     {
       if (num_cell_types > 1)
+      {
         throw std::runtime_error("Interior facet integrals with mixed "
                                  "topology aren't supported yet");
+      }
+
       const std::vector<int> c_offsets = a.coefficient_offsets();
       auto fn = a.kernel(IntegralType::interior_facet, i);
       assert(fn);

cpp/dolfinx/fem/assemble_vector_impl.h

@@ -1295,8 +1295,10 @@ void assemble_vector(
   std::shared_ptr<const mesh::Mesh<U>> mesh = L.mesh();
   assert(mesh);
   if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
+  {
     assemble_vector(b, L, mesh->geometry().dofmap(), mesh->geometry().x(),
                     constants, coefficients);
+  }
   else
   {
     auto x = mesh->geometry().x();