Remove gdim input to UFL elements, move value_shape to function space

cpp/dolfinx/fem/DirichletBC.h

@@ -335,7 +335,7 @@ class DirichletBC
   {
     assert(g);
     assert(V);
-    if (g->shape.size() != V->element()->value_shape().size())
+    if (g->shape.size() != V->value_shape().size())
     {
       throw std::runtime_error(
           "Rank mis-match between Constant and function space in DirichletBC");

cpp/dolfinx/fem/FiniteElement.cpp

@@ -123,7 +123,8 @@ _extract_sub_element(const FiniteElement<T>& finite_element,
 template <std::floating_point T>
 FiniteElement<T>::FiniteElement(const ufcx_finite_element& e)
     : _signature(e.signature), _space_dim(e.space_dimension),
-      _value_shape(e.value_shape, e.value_shape + e.value_rank),
+      _reference_value_shape(e.reference_value_shape,
+                             e.reference_value_shape + e.reference_value_rank),
       _bs(e.block_size), _is_mixed(e.element_type == ufcx_mixed_element)
 {
   const ufcx_shape _shape = e.cell_shape;
@@ -292,38 +293,20 @@ FiniteElement<T>::FiniteElement(const ufcx_finite_element& e)
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 FiniteElement<T>::FiniteElement(const basix::FiniteElement<T>& element,
-                                const std::vector<std::size_t>& value_shape)
-    : _value_shape(element.value_shape()), _is_mixed(false)
+                                const std::size_t block_size)
+    : _reference_value_shape(element.value_shape()), _bs(block_size),
+      _is_mixed(false)
 {
-  if (!_value_shape.empty() and !value_shape.empty())
-  {
-    throw std::runtime_error(
-        "Cannot specify value shape for non-scalar base element.");
-  }
-
-  // Set block size
-  if (!value_shape.empty())
-    _value_shape = value_shape;
-
-  // Set block size
-  if (!value_shape.empty())
-  {
-    _bs = std::accumulate(value_shape.begin(), value_shape.end(), 1,
-                          std::multiplies{});
-  }
-  else
-    _bs = 1;
-
   _space_dim = _bs * element.dim();
 
   // Create all sub-elements
   if (_bs > 1)
   {
     for (int i = 0; i < _bs; ++i)
     {
-      _sub_elements.push_back(std::make_shared<FiniteElement<T>>(
-          element, std::vector<std::size_t>{}));
+      _sub_elements.push_back(std::make_shared<FiniteElement<T>>(element, 1));
     }
+    _reference_value_shape = {block_size};
   }
 
   _element = std::make_unique<basix::FiniteElement<T>>(element);
@@ -388,17 +371,16 @@ int FiniteElement<T>::space_dimension() const noexcept
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-int FiniteElement<T>::value_size() const
+std::span<const std::size_t> FiniteElement<T>::reference_value_shape() const
 {
-  return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
-                         std::multiplies{});
+  return _reference_value_shape;
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
 int FiniteElement<T>::reference_value_size() const
 {
-  return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
-                         std::multiplies{});
+  return std::accumulate(_reference_value_shape.begin(),
+                         _reference_value_shape.end(), 1, std::multiplies{});
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
@@ -408,12 +390,6 @@ int FiniteElement<T>::block_size() const noexcept
 }
 //-----------------------------------------------------------------------------
 template <std::floating_point T>
-std::span<const std::size_t> FiniteElement<T>::value_shape() const noexcept
-{
-  return _value_shape;
-}
-//-----------------------------------------------------------------------------
-template <std::floating_point T>
 void FiniteElement<T>::tabulate(std::span<T> values, std::span<const T> X,
                                 std::array<std::size_t, 2> shape,
                                 int order) const

cpp/dolfinx/fem/FiniteElement.h

@@ -52,13 +52,9 @@ class FiniteElement
 
   /// @brief Create finite element from a Basix finite element.
   /// @param[in] element Basix finite element
-  /// @param[in] value_shape Value shape for 'blocked' elements, e.g.
-  /// vector-valued Lagrange elements where each component for the
-  /// vector field is a Lagrange element. For example, a vector-valued
-  /// element in 3D will have `value_shape` equal to `{3}`, and for a
-  /// second-order tensor element in 2D `value_shape` equal to `{2, 2}`.
+  /// @param[in] block_size The block size for the element
   FiniteElement(const basix::FiniteElement<geometry_type>& element,
-                const std::vector<std::size_t>& value_shape);
+                const std::size_t block_size);
 
   /// Copy constructor
   FiniteElement(const FiniteElement& element) = delete;
@@ -110,22 +106,13 @@ class FiniteElement
   /// @return Block size of the finite element space
   int block_size() const noexcept;
 
-  /// The value size, e.g. 1 for a scalar function, 2 for a 2D vector, 9
-  /// for a second-order tensor in 3D.
-  /// @note The return value of this function is equivalent to
-  /// `std::accumulate(value_shape().begin(), value_shape().end(), 1,
-  /// std::multiplies{})`.
-  /// @return The value size
-  int value_size() const;
-
   /// The value size, e.g. 1 for a scalar function, 2 for a 2D vector, 9
   /// for a second-order tensor in 3D, for the reference element
   /// @return The value size for the reference element
   int reference_value_size() const;
 
-  /// Shape of the value space. The rank is the size of the
-  /// `value_shape`.
-  std::span<const std::size_t> value_shape() const noexcept;
+  /// The reference value shape
+  std::span<const std::size_t> reference_value_shape() const;
 
   /// @brief Evaluate derivatives of the basis functions up to given order
   /// at points in the reference cell.
@@ -672,7 +659,7 @@ class FiniteElement
       _sub_elements;
 
   // Dimension of each value space
-  std::vector<std::size_t> _value_shape;
+  std::vector<std::size_t> _reference_value_shape;
 
   // Block size for BlockedElements. This gives the
   // number of DOFs co-located at each dof 'point'.

cpp/dolfinx/fem/Function.h

@@ -218,8 +218,7 @@ class Function
 
     const auto [fx, fshape] = f(_x);
     assert(fshape.size() <= 2);
-    if (int vs = _function_space->element()->value_size();
-        vs == 1 and fshape.size() == 1)
+    if (int vs = _function_space->value_size(); vs == 1 and fshape.size() == 1)
     {
       // Check for scalar-valued functions
       if (fshape.front() != x.size() / 3)
@@ -292,7 +291,7 @@ class Function
     if (e.argument_function_space())
       throw std::runtime_error("Cannot interpolate Expression with Argument");
 
-    if (value_size != _function_space->element()->value_size())
+    if (value_size != _function_space->value_size())
     {
       throw std::runtime_error(
           "Function value size not equal to Expression value size");
@@ -428,7 +427,7 @@ class Function
     const int bs_element = element->block_size();
     const std::size_t reference_value_size
         = element->reference_value_size() / bs_element;
-    const std::size_t value_size = element->value_size() / bs_element;
+    const std::size_t value_size = _function_space->value_size() / bs_element;
     const std::size_t space_dimension = element->space_dimension() / bs_element;
 
     // If the space has sub elements, concatenate the evaluations on the

cpp/dolfinx/fem/FunctionSpace.h

@@ -24,7 +24,6 @@
 
 namespace dolfinx::fem
 {
-
 /// @brief This class represents a finite element function space defined
 /// by a mesh, a finite element, and a local-to-global map of the
 /// degrees-of-freedom.
@@ -41,11 +40,14 @@ class FunctionSpace
   /// @param[in] mesh Mesh that the space is defined on.
   /// @param[in] element Finite element for the space.
   /// @param[in] dofmap Degree-of-freedom map for the space.
+  /// @param[in] value_shape The shape of the value space on the physical cell
   FunctionSpace(std::shared_ptr<const mesh::Mesh<geometry_type>> mesh,
                 std::shared_ptr<const FiniteElement<geometry_type>> element,
-                std::shared_ptr<const DofMap> dofmap)
+                std::shared_ptr<const DofMap> dofmap,
+                std::vector<std::size_t> value_shape)
       : _mesh(mesh), _element(element), _dofmap(dofmap),
-        _id(boost::uuids::random_generator()()), _root_space_id(_id)
+        _id(boost::uuids::random_generator()()), _root_space_id(_id),
+        _value_shape(value_shape)
   {
     // Do nothing
   }
@@ -91,7 +93,10 @@ class FunctionSpace
         = std::make_shared<DofMap>(_dofmap->extract_sub_dofmap(component));
 
     // Create new sub space
-    FunctionSpace sub_space(_mesh, element, dofmap);
+    FunctionSpace sub_space(_mesh, element, dofmap,
+                            compute_value_shape(element,
+                                                _mesh->topology()->dim(),
+                                                _mesh->geometry().dim()));
 
     // Set root space id and component w.r.t. root
     sub_space._root_space_id = _root_space_id;
@@ -150,8 +155,11 @@ class FunctionSpace
     auto collapsed_dofmap
         = std::make_shared<DofMap>(std::move(_collapsed_dofmap));
 
-    return {FunctionSpace(_mesh, _element, collapsed_dofmap),
-            std::move(collapsed_dofs)};
+    return {
+        FunctionSpace(_mesh, _element, collapsed_dofmap,
+                      compute_value_shape(_element, _mesh->topology()->dim(),
+                                          _mesh->geometry().dim())),
+        std::move(collapsed_dofs)};
   }
 
   /// @brief Get the component with respect to the root superspace.
@@ -316,6 +324,23 @@ class FunctionSpace
   /// The dofmap
   std::shared_ptr<const DofMap> dofmap() const { return _dofmap; }
 
+  /// The shape of the value space
+  std::span<const std::size_t> value_shape() const noexcept
+  {
+    return _value_shape;
+  }
+
+  /// The value size, e.g. 1 for a scalar-valued function, 2 for a 2D vector, 9
+  /// for a second-order tensor in 3D.
+  /// @note The return value of this function is equivalent to
+  /// `std::accumulate(value_shape().begin(), value_shape().end(), 1,
+  /// std::multiplies{})`.
+  int value_size() const
+  {
+    return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
+                           std::multiplies{});
+  }
+
 private:
   // The mesh
   std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
@@ -332,6 +357,8 @@ class FunctionSpace
   // Unique identifier for the space and for its root space
   boost::uuids::uuid _id;
   boost::uuids::uuid _root_space_id;
+
+  std::vector<std::size_t> _value_shape;
 };
 
 /// Extract FunctionSpaces for (0) rows blocks and (1) columns blocks
@@ -393,9 +420,41 @@ common_function_spaces(
   return {spaces0, spaces1};
 }
 
+/// @brief Compute the physical value shape of an element for a mesh
+/// @param[in] element The element
+/// @param[in] tdim Topological dimension
+/// @param[in] gdim Geometric dimension
+/// @return Physical valus shape
+template <std::floating_point T>
+std::vector<std::size_t> compute_value_shape(
+    std::shared_ptr<const dolfinx::fem::FiniteElement<T>> element,
+    std::size_t tdim, std::size_t gdim)
+{
+  auto rvs = element->reference_value_shape();
+  std::vector<std::size_t> value_shape(rvs.size());
+  if (element->block_size() > 1)
+  {
+    for (std::size_t i = 0; i < rvs.size(); ++i)
+    {
+      value_shape[i] = rvs[i];
+    }
+  }
+  else
+  {
+    for (std::size_t i = 0; i < rvs.size(); ++i)
+    {
+      if (rvs[i] == tdim)
+        value_shape[i] = gdim;
+      else
+        value_shape[i] = rvs[i];
+    }
+  }
+  return value_shape;
+}
+
 /// Type deduction
-template <typename U, typename V, typename W>
-FunctionSpace(U mesh, V element, W dofmap)
+template <typename U, typename V, typename W, typename X>
+FunctionSpace(U mesh, V element, W dofmap, X value_shape)
     -> FunctionSpace<typename std::remove_cvref<
         typename U::element_type>::type::geometry_type::value_type>;
 

cpp/dolfinx/fem/discreteoperators.h

@@ -197,8 +197,8 @@ void interpolation_matrix(const FunctionSpace<U>& V0,
   const std::size_t space_dim1 = e1->space_dimension();
   const std::size_t dim0 = space_dim0 / bs0;
   const std::size_t value_size_ref0 = e0->reference_value_size() / bs0;
-  const std::size_t value_size0 = e0->value_size() / bs0;
-  const std::size_t value_size1 = e1->value_size() / bs1;
+  const std::size_t value_size0 = V0.value_size() / bs0;
+  const std::size_t value_size1 = V1.value_size() / bs1;
 
   // Get geometry data
   const CoordinateElement<U>& cmap = mesh->geometry().cmap();
@@ -271,12 +271,12 @@ void interpolation_matrix(const FunctionSpace<U>& V0,
 
   // Basis values of Lagrange space unrolled for block size
   // (num_quadrature_points, Lagrange dof, value_size)
-  std::vector<U> basis_values_b(Xshape[0] * bs0 * dim0 * e1->value_size());
+  std::vector<U> basis_values_b(Xshape[0] * bs0 * dim0 * V1.value_size());
   mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0,
-                         e1->value_size());
-  std::vector<U> mapped_values_b(Xshape[0] * bs0 * dim0 * e1->value_size());
+                         V1.value_size());
+  std::vector<U> mapped_values_b(Xshape[0] * bs0 * dim0 * V1.value_size());
   mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0,
-                          e1->value_size());
+                          V1.value_size());
 
   auto pull_back_fn1
       = e1->basix_element().template map_fn<u_t, U_t, K_t, J_t>();
@@ -391,7 +391,7 @@ void interpolation_matrix(const FunctionSpace<U>& V0,
     {
       MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
           T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
-          A(Ab.data(), Xshape[0], e1->value_size(), space_dim0);
+          A(Ab.data(), Xshape[0], V1.value_size(), space_dim0);
       for (std::size_t i = 0; i < mapped_values.extent(0); ++i)
         for (std::size_t j = 0; j < mapped_values.extent(1); ++j)
           for (std::size_t k = 0; k < mapped_values.extent(2); ++k)