merge latest changes of main into sclaus branch

Author: Claus Susanne

cpp/dolfinx/fem/CMakeLists.txt

@@ -10,6 +10,7 @@ set(HEADERS_fem
     ${CMAKE_CURRENT_SOURCE_DIR}/Function.h
     ${CMAKE_CURRENT_SOURCE_DIR}/FunctionSpace.h
     ${CMAKE_CURRENT_SOURCE_DIR}/assembler.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/assemble_expression_impl.h
     ${CMAKE_CURRENT_SOURCE_DIR}/assemble_matrix_impl.h
     ${CMAKE_CURRENT_SOURCE_DIR}/assemble_scalar_impl.h
     ${CMAKE_CURRENT_SOURCE_DIR}/assemble_vector_impl.h

cpp/dolfinx/fem/Expression.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2020-2021 Jack S. Hale and Michal Habera.
+// Copyright (C) 2020-2025 Jack S. Hale, Michal Habera and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -8,32 +8,34 @@
 
 #include "Constant.h"
 #include "Function.h"
-#include "pack.h"
-#include "utils.h"
 #include <algorithm>
 #include <array>
 #include <concepts>
 #include <dolfinx/common/types.h>
 #include <dolfinx/mesh/Mesh.h>
 #include <functional>
+#include <memory>
 #include <span>
 #include <utility>
 #include <vector>
 
 namespace dolfinx::fem
 {
-/// @brief Represents a mathematical expression evaluated at a
-/// pre-defined set of points on the reference cell.
+template <dolfinx::scalar T, std::floating_point U>
+class Function;
+
+/// @brief An Expression represents a mathematical expression evaluated
+/// at a pre-defined points on a reference cell.
 ///
-//// This class closely follows the concept of a UFC Expression.
+/// An Expression can be evaluated using ::tabulate_expression.
 ///
-/// The functionality can be used to evaluate a gradient of a Function
-/// at quadrature points in all cells. This evaluated gradient can then
-/// be used as input in to a non-FEniCS function that calculates a
+/// An example of Expression use is to evaluate the gradient of a
+/// Function at quadrature points in cells. The evaluated gradient can
+/// then be used as input in to a non-FEniCS function that evaluates a
 /// material constitutive model.
 ///
-/// @tparam T The scalar type
-/// @tparam U The mesh geometry scalar type
+/// @tparam T The scalar type.
+/// @tparam U The mesh geometry scalar type.
 template <dolfinx::scalar T,
           std::floating_point U = dolfinx::scalar_value_type_t<T>>
 class Expression
@@ -50,18 +52,21 @@ class Expression
 
   /// @brief Create an Expression.
   ///
-  /// @note Users should prefer the @ref create_expression factory
-  /// functions.
+  /// @note Users should prefer the fem::create_expression factory
+  /// functions for creating an Expression.
   ///
   /// @param[in] coefficients Coefficients in the Expression.
-  /// @param[in] constants Constants in the Expression
+  /// @param[in] constants Constants in the Expression.
   /// @param[in] X Points on the reference cell, `shape=(number of
   /// points, tdim)` and storage is row-major.
   /// @param[in] Xshape Shape of `X`.
   /// @param[in] fn Function for tabulating the Expression.
   /// @param[in] value_shape Shape of Expression evaluated at single
   /// point.
-  /// @param[in] argument_function_space Function space for Argument.
+  /// @param[in] argument_space Function space for Argument. Used to
+  /// computed a 1-form expression, e.g. can be used to create a matrix
+  /// that when applied to a degree-of-freedom vector gives the
+  /// expression values at the evaluation points.
   Expression(
       const std::vector<std::shared_ptr<
           const Function<scalar_type, geometry_type>>>& coefficients,
@@ -72,13 +77,11 @@ class Expression
                          const geometry_type*, const int*, const uint8_t*, void*)>
           fn,
       const std::vector<std::size_t>& value_shape,
-      std::shared_ptr<const FunctionSpace<geometry_type>>
-          argument_function_space
+      std::shared_ptr<const FunctionSpace<geometry_type>> argument_space
       = nullptr)
       : _coefficients(coefficients), _constants(constants),
         _x_ref(std::vector<geometry_type>(X.begin(), X.end()), Xshape), _fn(fn),
-        _value_shape(value_shape),
-        _argument_function_space(argument_function_space)
+        _value_shape(value_shape), _argument_space(argument_space)
   {
     for (auto& c : _coefficients)
     {
@@ -97,28 +100,25 @@ class Expression
   /// Destructor
   virtual ~Expression() = default;
 
-  /// @brief Get argument function space.
-  /// @return The argument function space, nullptr if there is no
-  /// argument.
-  std::shared_ptr<const FunctionSpace<geometry_type>>
-  argument_function_space() const
+  /// @brief Argument function space.
+  /// @return Argument function space, nullptr if there is no argument.
+  std::shared_ptr<const FunctionSpace<geometry_type>> argument_space() const
   {
-    return _argument_function_space;
+    return _argument_space;
   };
 
-  /// @brief Get coefficients.
-  /// @return Vector of attached coefficients.
+  /// @brief Expression coefficients.
+  /// @return List of attached coefficients.
   const std::vector<
       std::shared_ptr<const Function<scalar_type, geometry_type>>>&
   coefficients() const
   {
     return _coefficients;
   }
 
-  /// @brief Get constants.
-  /// @return Vector of attached constants with their names. Names are
-  /// used to set constants in user's c++ code. Index in the vector is
-  /// the position of the constant in the original (nonsimplified) form.
+  /// @brief Expression constants.
+  ///
+  /// @return List of attached constants.
   const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
   constants() const
   {
@@ -127,7 +127,7 @@ class Expression
 
   /// @brief Offset for each coefficient expansion array on a cell.
   ///
-  /// Used to pack data for multiple coefficients in a flat array. The
+  /// Used to pack data for multiple coefficients into a flat array. The
   /// last entry is the size required to store all coefficients.
   /// @return The offsets.
   std::vector<int> coefficient_offsets() const
@@ -142,157 +142,33 @@ class Expression
     return n;
   }
 
-  /// @brief Evaluate Expression on cells or facets.
-  ///
-  /// @param[in] mesh Cells on which to evaluate the Expression.
-  /// @param[in] entities List of entities to evaluate the expression
-  /// on. This could be either a list of cells or a list of (cell, local
-  /// @param[out] values A 2D array to store the result. Caller is
-  /// responsible for correct sizing which should be `(num_cells,
-  /// num_points * value_size * num_all_argument_dofs columns)`.
-  /// facet index) tuples. Array is flattened per entity.
-  /// @param[in] vshape The shape of `values` (row-major storage).
-  void eval(const mesh::Mesh<geometry_type>& mesh,
-            std::span<const std::int32_t> entities,
-            std::span<scalar_type> values,
-            std::array<std::size_t, 2> vshape) const
-  {
-    std::size_t estride;
-    if (mesh.topology()->dim() == _x_ref.second[1])
-      estride = 1;
-    else if (mesh.topology()->dim() == _x_ref.second[1] + 1)
-      estride = 2;
-    else
-      throw std::runtime_error("Invalid dimension of evaluation points.");
-
-    // Prepare coefficients and constants
-    std::vector<T> coeffs((entities.size() / estride)
-                          * this->coefficient_offsets().back());
-    int cstride = this->coefficient_offsets().back();
-    {
-      std::vector<
-          std::reference_wrapper<const Function<scalar_type, geometry_type>>>
-          c;
-      std::ranges::transform(this->coefficients(), std::back_inserter(c),
-                             [](auto c) -> const Function<T, U>&
-                             { return *c; });
-      fem::pack_coefficients(c, this->coefficient_offsets(), entities, estride,
-                             std::span(coeffs));
-    }
-    std::vector<scalar_type> constant_data = fem::pack_constants(*this);
-
-    auto fn = this->get_tabulate_expression();
-
-    // Prepare cell geometry
-    auto x_dofmap = mesh.geometry().dofmap();
-
-    // Get geometry data
-    auto& cmap = mesh.geometry().cmap();
-
-    std::size_t num_dofs_g = cmap.dim();
-    auto x_g = mesh.geometry().x();
-
-    // Create data structures used in evaluation
-    std::vector<geometry_type> coord_dofs(3 * num_dofs_g);
-
-    int num_argument_dofs = 1;
-    std::span<const std::uint32_t> cell_info;
-    std::function<void(std::span<scalar_type>, std::span<const std::uint32_t>,
-                       std::int32_t, int)>
-        post_dof_transform
-        = [](std::span<scalar_type>, std::span<const std::uint32_t>,
-             std::int32_t, int)
-    {
-      // Do nothing
-    };
-
-    if (_argument_function_space)
-    {
-      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())
-      {
-        mesh.topology_mutable()->create_entity_permutations();
-        cell_info = std::span(mesh.topology()->get_cell_permutation_info());
-        post_dof_transform
-            = element->template dof_transformation_right_fn<scalar_type>(
-                doftransform::transpose);
-      }
-    }
-
-    // Create get entity index function
-    std::function<const std::int32_t*(std::span<const std::int32_t>,
-                                      std::size_t)>
-        get_entity_index
-        = []([[maybe_unused]] std::span<const std::int32_t> entities,
-             [[maybe_unused]] std::size_t idx) { return nullptr; };
-    if (estride == 2)
-    {
-      get_entity_index
-          = [](std::span<const std::int32_t> entities, std::size_t idx)
-      { return entities.data() + 2 * idx + 1; };
-    }
-
-    // Iterate over cells and 'assemble' into values
-    int size0 = _x_ref.second[0] * value_size();
-    std::vector<scalar_type> values_local(size0 * num_argument_dofs, 0);
-    for (std::size_t e = 0; e < entities.size() / estride; ++e)
-    {
-      std::int32_t entity = entities[e * estride];
-      auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
-          x_dofmap, entity, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
-      for (std::size_t i = 0; i < x_dofs.size(); ++i)
-      {
-        std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3,
-                    std::next(coord_dofs.begin(), 3 * i));
-      }
-
-      const scalar_type* coeff_cell = coeffs.data() + e * cstride;
-      const int* entity_index = get_entity_index(entities, e);
-
-      std::ranges::fill(values_local, 0);
-      _fn(values_local.data(), coeff_cell, constant_data.data(),
-          coord_dofs.data(), entity_index, nullptr, 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];
-    }
-  }
-
-  /// @brief Get function for tabulate_expression.
-  /// @return fn Function to tabulate expression.
+  /// @brief Function for tabulating the Expression.
   const std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
                            const geometry_type*, const int*, const uint8_t*, void*)>&
-  get_tabulate_expression() const
+  kernel() const
   {
     return _fn;
   }
 
-  /// @brief Get value size
-  /// @return The value size.
+  /// @brief Value size of the Expression result.
   int value_size() const
   {
     return std::reduce(_value_shape.begin(), _value_shape.end(), 1,
                        std::multiplies{});
   }
 
-  /// @brief Get value shape.
-  /// @return The value shape.
+  /// @brief Value shape of of Expression result (at a point),
   const std::vector<std::size_t>& value_shape() const { return _value_shape; }
 
-  /// @brief Evaluation points on the reference cell.
-  /// @return Evaluation points.
+  /// @brief Evaluation point coordinates on the reference cell.
   std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> X() const
   {
     return _x_ref;
   }
 
 private:
   // Function space for Argument
-  std::shared_ptr<const FunctionSpace<geometry_type>> _argument_function_space;
+  std::shared_ptr<const FunctionSpace<geometry_type>> _argument_space;
 
   // Coefficients associated with the Expression
   std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>>
@@ -306,11 +182,10 @@ class Expression
                      const geometry_type*, const int*, const uint8_t*, void*)>
       _fn;
 
-  // Shape of the evaluated expression
+  // Shape of the evaluated Expression
   std::vector<std::size_t> _value_shape;
 
-  // Evaluation points on reference cell. Synonymous with X in public
-  // interface.
+  // Evaluation points on reference cell
   std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _x_ref;
 };
 } // namespace dolfinx::fem

cpp/dolfinx/fem/Function.h

@@ -9,8 +9,10 @@
 #include "DofMap.h"
 #include "FiniteElement.h"
 #include "FunctionSpace.h"
+#include "assembler.h"
 #include "interpolate.h"
 #include <algorithm>
+#include <basix/mdspan.hpp>
 #include <concepts>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/common/types.h>
@@ -315,6 +317,8 @@ class Function
                    std::span<const std::int32_t> cells0,
                    std::span<const std::int32_t> cells1 = {})
   {
+    namespace md = MDSPAN_IMPL_STANDARD_NAMESPACE;
+
     // Extract mesh
     const mesh::Mesh<geometry_type>* mesh0 = nullptr;
     for (auto& c : e0.coefficients())
@@ -352,7 +356,7 @@ class Function
     // Check that Function and Expression spaces are compatible
     assert(_function_space->element());
     std::size_t value_size = e0.value_size();
-    if (e0.argument_function_space())
+    if (e0.argument_space())
       throw std::runtime_error("Cannot interpolate Expression with Argument.");
     if (value_size != _function_space->element()->value_size())
     {
@@ -391,7 +395,8 @@ class Function
         f(fdata.data(), num_cells, num_points, value_size);
 
     // Evaluate Expression at points
-    e0.eval(*mesh0, cells0, fdata, {num_cells, num_points * value_size});
+    tabulate_expression(std::span(fdata), e0, *mesh0,
+                        md::mdspan(cells0.data(), cells0.size()));
 
     // Reshape evaluated data to fit interpolate.
     // Expression returns matrix of shape (num_cells, num_points *