Vertex integrals

.github/workflows/ccpp.yml

@@ -97,9 +97,9 @@ jobs:
 
       - name: Install FEniCS Python components
         run: |
-          pip install --no-build-isolation git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }}
-          pip install --no-build-isolation git+https://github.com/fenics/basix.git@${{ env.basix_ref }}
-          pip install --no-build-isolation git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
 
       - name: Configure (C++)
         working-directory: cpp
@@ -189,9 +189,9 @@ jobs:
 
       - name: Install FEniCS Python components
         run: |
-          pip install --no-build-isolation git+https://github.com/FEniCS/ufl.git@${{ env.ufl_ref }}
-          pip install --no-build-isolation git+https://github.com/FEniCS/basix.git@${{ env.basix_ref }}
-          pip install --no-build-isolation git+https://github.com/FEniCS/ffcx.git@${{ env.ffcx_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
 
       - name: Configure (C++)
         working-directory: cpp
@@ -293,9 +293,9 @@ jobs:
       - name: Install FEniCS Python components
         run: |
           pip install -r python/build-requirements.txt
-          pip install --no-build-isolation git+https://github.com/FEniCS/ufl.git@${{ env.ufl_ref }}
-          pip install --no-build-isolation git+https://github.com/FEniCS/basix.git@${{ env.basix_ref }}
-          pip install --no-build-isolation git+https://github.com/FEniCS/ffcx.git@${{ env.ffcx_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
 
       - name: Configure C++
         run: |

.github/workflows/ci-spack.yml

@@ -86,9 +86,9 @@ jobs:
         run: |
           . ./spack-src/share/spack/setup-env.sh
           spack env activate .
-          pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }}
-          pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }}
-          pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }}
+          pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
 
       - name: Configure and build C++
         run: |

.github/workflows/fenicsx-refs.env

@@ -1,4 +1,6 @@
+basix_repository=FEniCS/basix
 basix_ref=main
+ufl_repository=FEniCS/ufl
 ufl_ref=main
+ffcx_repository=FEniCS/ffcx
 ffcx_ref=main
-

.github/workflows/macos.yml

@@ -73,9 +73,9 @@ jobs:
 
       - name: Install FEniCSx dependencies
         run: |
-          python -m pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }}
-          python -m pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }}
-          python -m pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }}
+          python -m pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          python -m pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          python -m pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
 
       - name: Configure (C++)
         working-directory: cpp

.github/workflows/sonarcloud.yml

@@ -62,9 +62,9 @@ jobs:
           echo "$HOME/.sonar/build-wrapper-linux-x86" >> $GITHUB_PATH
       - name: Install FEniCS Python components
         run: |
-          python -m pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }}
-          python -m pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }}
-          python -m pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }}
+          pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }}
+          pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }}
+          pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }}
       - name: Run build-wrapper
         run: |
           mkdir build

cpp/dolfinx/fem/CMakeLists.txt

@@ -18,6 +18,7 @@ set(HEADERS_fem
     ${CMAKE_CURRENT_SOURCE_DIR}/dofmapbuilder.h
     ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_fem.h
     ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/kernel.h
     ${CMAKE_CURRENT_SOURCE_DIR}/petsc.h
     ${CMAKE_CURRENT_SOURCE_DIR}/pack.h
     ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.h

cpp/dolfinx/fem/assemble_scalar_impl.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2019-2025 Garth N. Wells
+// Copyright (C) 2019-2025 Garth N. Wells and Paul T. Kühner
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -150,6 +150,43 @@ T assemble_interior_facets(
   return value;
 }
 
+/// Assemble functional over vertices
+template <dolfinx::scalar T>
+T assemble_vertices(mdspan2_t x_dofmap,
+                    md::mdspan<const scalar_value_t<T>,
+                               md::extents<std::size_t, md::dynamic_extent, 3>>
+                        x,
+                    md::mdspan<const std::int32_t,
+                               md::extents<std::size_t, md::dynamic_extent, 2>>
+                        vertices,
+                    FEkernel<T> auto fn, std::span<const T> constants,
+                    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs)
+{
+  T value(0);
+  if (vertices.empty())
+    return value;
+
+  // Create data structures used in assembly
+  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+
+  // Iterate over all cells
+  for (std::size_t index = 0; index < vertices.extent(0); ++index)
+  {
+    std::int32_t cell = vertices(index, 0);
+    std::int32_t local_vertex_index = vertices(index, 1);
+
+    // Get cell coordinates/geometry
+    auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
+    for (std::size_t i = 0; i < x_dofs.size(); ++i)
+      std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
+
+    fn(&value, &coeffs(index, 0), constants.data(), cdofs.data(),
+       &local_vertex_index, nullptr, nullptr);
+  }
+
+  return value;
+}
+
 /// Assemble functional into an scalar with provided mesh geometry.
 template <dolfinx::scalar T, std::floating_point U>
 T assemble_scalar(
@@ -198,14 +235,19 @@ T assemble_scalar(
         = coefficients.at({IntegralType::exterior_facet, i});
 
     std::span facets = M.domain(IntegralType::exterior_facet, i, 0);
+
+    constexpr std::size_t num_adjacent_cells = 1;
+    // Two values per each adj. cell (cell index and local facet index).
+    constexpr std::size_t shape1 = 2 * num_adjacent_cells;
+
     assert((facets.size() / 2) * cstride == coeffs.size());
     value += impl::assemble_exterior_facets(
         x_dofmap, x,
         md::mdspan<const std::int32_t,
                    md::extents<std::size_t, md::dynamic_extent, 2>>(
-            facets.data(), facets.size() / 2, 2),
-        fn, constants, md::mdspan(coeffs.data(), facets.size() / 2, cstride),
-        perms);
+            facets.data(), facets.size() / shape1, 2),
+        fn, constants,
+        md::mdspan(coeffs.data(), facets.size() / shape1, cstride), perms);
   }
 
   for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i)
@@ -215,19 +257,48 @@ T assemble_scalar(
     auto& [coeffs, cstride]
         = coefficients.at({IntegralType::interior_facet, i});
     std::span facets = M.domain(IntegralType::interior_facet, i, 0);
-    assert((facets.size() / 4) * 2 * cstride == coeffs.size());
+
+    constexpr std::size_t num_adjacent_cells = 2;
+    // Two values per each adj. cell (cell index and local facet index).
+    constexpr std::size_t shape1 = 2 * num_adjacent_cells;
+
+    assert((facets.size() / shape1) * 2 * cstride == coeffs.size());
     value += impl::assemble_interior_facets(
         x_dofmap, x,
         md::mdspan<const std::int32_t,
                    md::extents<std::size_t, md::dynamic_extent, 2, 2>>(
-            facets.data(), facets.size() / 4, 2, 2),
+            facets.data(), facets.size() / shape1, 2, 2),
         fn, constants,
         md::mdspan<const T, md::extents<std::size_t, md::dynamic_extent, 2,
                                         md::dynamic_extent>>(
-            coeffs.data(), facets.size() / 4, 2, cstride),
+            coeffs.data(), facets.size() / shape1, 2, cstride),
         perms);
   }
 
+  for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i)
+  {
+    auto fn = M.kernel(IntegralType::vertex, i, 0);
+    assert(fn);
+
+    auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
+
+    std::span<const std::int32_t> vertices
+        = M.domain(IntegralType::vertex, i, 0);
+    assert(vertices.size() * cstride == coeffs.size());
+
+    constexpr std::size_t num_adjacent_cells = 1;
+    // Two values per adj. cell (cell index and local vertex index).
+    constexpr std::size_t shape1 = 2 * num_adjacent_cells;
+
+    value += impl::assemble_vertices(
+        x_dofmap, x,
+        md::mdspan<const std::int32_t,
+                   md::extents<std::size_t, md::dynamic_extent, 2>>(
+            vertices.data(), vertices.size() / shape1, shape1),
+        fn, constants,
+        md::mdspan(coeffs.data(), vertices.size() / shape1, cstride));
+  }
+
   return value;
 }
 

cpp/dolfinx/fem/assemble_vector_impl.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2018-2021 Garth N. Wells
+// Copyright (C) 2018-2025 Garth N. Wells and Paul T. Kühner
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -941,6 +941,94 @@ void assemble_interior_facets(
   }
 }
 
+/// @brief Execute kernel over a set of vertices and accumulate result in
+/// vector.
+///
+/// @tparam T  Scalar type
+/// @tparam _bs Block size of the form test function dof map. If less
+/// than zero the block size is determined at runtime. If `_bs` is
+/// positive the block size is used as a compile-time constant, which
+/// has performance benefits.
+/// @param[in] P0 Function that applies transformation `P0.b` in-place
+/// to `b` to transform test degrees-of-freedom.
+/// @param[in,out] b Array to accumulate into.
+/// @param[in] x_dofmap Dofmap for the mesh geometry.
+/// @param[in] x Mesh geometry (coordinates).
+/// @param[in] vertices Vertex indices `(vertices.size(), 2)` - first entry
+/// holds the index of the cell adjacent to the vertex, and the second
+/// stores the local index of the vertex within the cell.
+/// @param[in] dofmap Test function (row) degree-of-freedom data holding
+/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices.
+/// @param[in] kernel Kernel function to execute over each cell.
+/// @param[in] constants Constant coefficient data in the kernel.
+/// @param[in] coeffs Coefficient data in the kernel. It has shape
+/// `(vertices.size(), num_cell_coeffs)`. `coeffs(i, j)` is the `j`th
+/// coefficient for cell `i`.
+/// @param[in] cell_info0 Cell permutation information for the test
+/// function mesh.
+template <dolfinx::scalar T, int _bs = -1>
+void assemble_vertices(
+    fem::DofTransformKernel<T> auto P0, std::span<T> b, mdspan2_t x_dofmap,
+    md::mdspan<const scalar_value_t<T>,
+               md::extents<std::size_t, md::dynamic_extent, 3>>
+        x,
+    md::mdspan<const std::int32_t,
+               md::extents<std::size_t, md::dynamic_extent, 2>>
+        vertices,
+    std::tuple<mdspan2_t, int,
+               md::mdspan<const std::int32_t,
+                          md::extents<std::size_t, md::dynamic_extent, 2>>>
+        dofmap,
+    FEkernel<T> auto kernel, std::span<const T> constants,
+    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
+    std::span<const std::uint32_t> cell_info0)
+{
+  if (vertices.empty())
+    return;
+
+  const auto [dmap, bs, vertices0] = dofmap;
+  assert(_bs < 0 or _bs == bs);
+
+  // Create data structures used in assembly
+  std::vector<scalar_value_t<T>> cdofs(3 * x_dofmap.extent(1));
+  std::vector<T> be(bs * dmap.extent(1));
+
+  // Iterate over active vertices
+  for (std::size_t index = 0; index < vertices.extent(0); ++index)
+  {
+    // Integration domain cell, local index, and test function cell
+    std::int32_t cell = vertices(index, 0);
+    std::int32_t local_index = vertices(index, 1);
+    std::int32_t c0 = vertices0(index, 0);
+
+    // Get cell coordinates/geometry
+    auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent);
+    for (std::size_t i = 0; i < x_dofs.size(); ++i)
+      std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i));
+
+    // Tabulate vector for vertex
+    std::ranges::fill(be, 0);
+    kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(),
+           &local_index, nullptr, nullptr);
+    P0(be, cell_info0, c0, 1);
+
+    // Scatter vertex vector to 'global' vector array
+    auto dofs = md::submdspan(dmap, c0, md::full_extent);
+    if constexpr (_bs > 0)
+    {
+      for (std::size_t i = 0; i < dofs.size(); ++i)
+        for (int k = 0; k < _bs; ++k)
+          b[_bs * dofs[i] + k] += be[_bs * i + k];
+    }
+    else
+    {
+      for (std::size_t i = 0; i < dofs.size(); ++i)
+        for (int k = 0; k < bs; ++k)
+          b[bs * dofs[i] + k] += be[bs * i + k];
+    }
+  }
+}
+
 /// Modify RHS vector to account for boundary condition such that:
 ///
 /// b <- b - alpha * A.(x_bc - x0)
@@ -1381,6 +1469,32 @@ void assemble_vector(
             cell_info0, perms);
       }
     }
+
+    for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i)
+    {
+      auto fn = L.kernel(IntegralType::vertex, i, 0);
+      assert(fn);
+
+      std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx);
+      std::span vertices0
+          = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx);
+
+      auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i});
+
+      assert(vertices.size() * cstride == coeffs.size());
+
+      impl::assemble_vertices<T>(
+          P0, b, x_dofmap, x,
+          md::mdspan<const std::int32_t,
+                     md::extents<std::size_t, md::dynamic_extent, 2>>(
+              vertices.data(), vertices.size() / 2, 2),
+          {dofs, bs,
+           md::mdspan<const std::int32_t,
+                      md::extents<std::size_t, md::dynamic_extent, 2>>(
+               vertices0.data(), vertices0.size() / 2, 2)},
+          fn, constants,
+          md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0);
+    }
   }
 }