Encourage proper management of subspaces in user code (#2916)

* return new object not share pointer * update demo * encourage users to hold the subspace * add note to documentation * remove dot * return a FunctionSpace instead of shared pointer * Fix template type * Formatting fix * Update Pythin docs * Simplifications * Template type updates * Doc fix --------- Co-authored-by: Garth N. Wells <[email protected]>
FEniCS · Dec 1, 2023 · 30d09ec · 30d09ec
1 parent 8eefd77
commit 30d09ec
Show file tree

Hide file tree

Showing 8 changed files with 110 additions and 81 deletions.
diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
@@ -205,7 +205,7 @@ int main(int argc, char* argv[])
           }
           return marker;
         });
-    auto bcs = std::vector{
+    std::vector bcs = {
         std::make_shared<const fem::DirichletBC<T>>(std::vector<T>{0, 0, 0},
                                                     bdofs_left, V),
         std::make_shared<const fem::DirichletBC<T>>(u_rotation, bdofs_right)};

diff --git a/cpp/dolfinx/fem/Constant.h b/cpp/dolfinx/fem/Constant.h
@@ -22,27 +22,30 @@ template <dolfinx::scalar T>
 class Constant
 {
 public:
+  /// Field type
+  using value_type = T;
+
   /// @brief Create a rank-0 (scalar-valued) constant
   /// @param[in] c Value of the constant
-  explicit Constant(T c) : value({c}) {}
+  explicit Constant(value_type c) : value({c}) {}
 
   /// @brief Create a rank-1 (vector-valued) constant
   /// @param[in] c Value of the constant
-  explicit Constant(std::span<const T> c)
+  explicit Constant(std::span<const value_type> c)
       : Constant(c, std::vector<std::size_t>{c.size()})
   {
   }
 
   /// @brief Create a rank-d constant
   /// @param[in] c Value of the Constant (row-majors storage)
   /// @param[in] shape Shape of the Constant
-  Constant(std::span<const T> c, std::span<const std::size_t> shape)
+  Constant(std::span<const value_type> c, std::span<const std::size_t> shape)
       : value(c.begin(), c.end()), shape(shape.begin(), shape.end())
   {
   }
 
   /// Values, stored as a row-major flattened array
-  std::vector<T> value;
+  std::vector<value_type> value;
 
   /// Shape
   std::vector<std::size_t> shape;

diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h
@@ -26,6 +26,9 @@ template <std::floating_point T>
 class FiniteElement
 {
 public:
+  /// Geometry type of the Mesh that the FunctionSpace is defined on.
+  using geometry_type = T;
+
   /// @brief Create finite element from UFC finite element.
   /// @param[in] e UFC finite element.
   explicit FiniteElement(const ufcx_finite_element& e);
@@ -37,7 +40,7 @@ class FiniteElement
   /// 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}`.
-  FiniteElement(const basix::FiniteElement<T>& element,
+  FiniteElement(const basix::FiniteElement<geometry_type>& element,
                 const std::vector<std::size_t>& value_shape);
 
   /// Copy constructor
@@ -118,7 +121,8 @@ class FiniteElement
   /// @param[in] shape The shape of `X`
   /// @param[in] order The number of derivatives (up to and including
   /// this order) to tabulate for
-  void tabulate(std::span<T> values, std::span<const T> X,
+  void tabulate(std::span<geometry_type> values,
+                std::span<const geometry_type> X,
                 std::array<std::size_t, 2> shape, int order) const;
 
   /// Evaluate all derivatives of the basis functions up to given order
@@ -130,8 +134,8 @@ class FiniteElement
   /// @param[in] order The number of derivatives (up to and including
   /// this order) to tabulate for
   /// @return Basis function values and array shape (row-major storage)
-  std::pair<std::vector<T>, std::array<std::size_t, 4>>
-  tabulate(std::span<const T> X, std::array<std::size_t, 2> shape,
+  std::pair<std::vector<geometry_type>, std::array<std::size_t, 4>>
+  tabulate(std::span<const geometry_type> X, std::array<std::size_t, 2> shape,
            int order) const;
 
   /// @brief Number of sub elements (for a mixed or blocked element)
@@ -145,15 +149,15 @@ class FiniteElement
   bool is_mixed() const noexcept;
 
   /// Subelements (if any)
-  const std::vector<std::shared_ptr<const FiniteElement<T>>>&
+  const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
   sub_elements() const noexcept;
 
   /// Extract sub finite element for component
-  std::shared_ptr<const FiniteElement<T>>
+  std::shared_ptr<const FiniteElement<geometry_type>>
   extract_sub_element(const std::vector<int>& component) const;
 
   /// Return underlying basix element (if it exists)
-  const basix::FiniteElement<T>& basix_element() const;
+  const basix::FiniteElement<geometry_type>& basix_element() const;
 
   /// Get the map type used by the element
   basix::maps::type map_type() const;
@@ -182,7 +186,7 @@ class FiniteElement
   /// @return Interpolation point coordinates on the reference cell,
   /// returning the (0) coordinates data (row-major) storage and (1) the
   /// shape `(num_points, tdim)`.
-  std::pair<std::vector<T>, std::array<std::size_t, 2>>
+  std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
   interpolation_points() const;
 
   /// Interpolation operator (matrix) `Pi` that maps a function
@@ -194,7 +198,7 @@ class FiniteElement
   /// @return The interpolation operator `Pi`, returning the data for
   /// `Pi` (row-major storage) and the shape `(num_dofs, num_points *
   /// value_size)`
-  std::pair<std::vector<T>, std::array<std::size_t, 2>>
+  std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
   interpolation_operator() const;
 
   /// @brief Create a matrix that maps degrees of freedom from one
@@ -209,7 +213,7 @@ class FiniteElement
   /// @pre The two elements must use the same mapping between the
   /// reference and physical cells
   /// @note Does not support mixed elements
-  std::pair<std::vector<T>, std::array<std::size_t, 2>>
+  std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
   create_interpolation_operator(const FiniteElement& from) const;
 
   /// @brief Check if DOF transformations are needed for this element.
@@ -671,7 +675,8 @@ class FiniteElement
   int _space_dim;
 
   // List of sub-elements (if any)
-  std::vector<std::shared_ptr<const FiniteElement<T>>> _sub_elements;
+  std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>
+      _sub_elements;
 
   // Dimension of each value space
   std::vector<std::size_t> _value_shape;
@@ -685,10 +690,10 @@ class FiniteElement
   bool _needs_dof_transformations;
 
   // Basix Element (nullptr for mixed elements)
-  std::unique_ptr<basix::FiniteElement<T>> _element;
+  std::unique_ptr<basix::FiniteElement<geometry_type>> _element;
 
   // Quadrature points of a quadrature element (0 dimensional array for
   // all elements except quadrature elements)
-  std::pair<std::vector<T>, std::array<std::size_t, 2>> _points;
+  std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _points;
 };
 } // namespace dolfinx::fem
diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h
@@ -100,12 +100,13 @@ class Function
   // Assignment
   Function& operator=(const Function& v) = delete;
 
-  /// @brief Extract sub-function (a view into the Function).
+  /// @brief Extract a sub-function (a view into the Function).
   /// @param[in] i Index of subfunction
   /// @return The sub-function
   Function sub(int i) const
   {
-    auto sub_space = _function_space->sub({i});
+    auto sub_space = std::make_shared<FunctionSpace<geometry_type>>(
+        _function_space->sub({i}));
     assert(sub_space);
     return Function(sub_space, _x);
   }

diff --git a/cpp/dolfinx/fem/FunctionSpace.h b/cpp/dolfinx/fem/FunctionSpace.h
@@ -28,18 +28,21 @@ 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.
-/// @tparam T The floating point (real) type of the mesh geometry and the
-/// finite element basis.
+/// @tparam T The floating point (real) type of the mesh geometry and
+/// the finite element basis.
 template <std::floating_point T>
 class FunctionSpace
 {
 public:
+  /// Geometry type of the Mesh that the FunctionSpace is defined on.
+  using geometry_type = T;
+
   /// @brief Create function space for given mesh, element and dofmap.
   /// @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.
-  FunctionSpace(std::shared_ptr<const mesh::Mesh<T>> mesh,
-                std::shared_ptr<const FiniteElement<T>> element,
+  FunctionSpace(std::shared_ptr<const mesh::Mesh<geometry_type>> mesh,
+                std::shared_ptr<const FiniteElement<geometry_type>> element,
                 std::shared_ptr<const DofMap> dofmap)
       : _mesh(mesh), _element(element), _dofmap(dofmap),
         _id(boost::uuids::random_generator()()), _root_space_id(_id)
@@ -62,21 +65,23 @@ class FunctionSpace
   /// Move assignment operator
   FunctionSpace& operator=(FunctionSpace&& V) = default;
 
-  /// Extract subspace for component
-  /// @param[in] component The subspace component
-  /// @return The subspace
-  std::shared_ptr<FunctionSpace<T>> sub(const std::vector<int>& component) const
+  ///  @brief Create a subspace (view) for a specific component.
+  ///
+  /// @note If the subspace is re-used, for performance reasons the
+  /// returned subspace should be stored by the caller to avoid repeated
+  /// re-computation of the subspace.
+  ///
+  /// @param[in] component Subspace component.
+  /// @return A subspace.
+  FunctionSpace sub(const std::vector<int>& component) const
   {
     assert(_mesh);
     assert(_element);
     assert(_dofmap);
 
-    // Check if sub space is already in the cache and not expired
-    if (auto it = _subspaces.find(component); it != _subspaces.end())
-    {
-      if (auto s = it->second.lock())
-        return s;
-    }
+    // Check that component is valid
+    if (component.empty())
+      throw std::runtime_error("Component must be non-empty");
 
     // Extract sub-element
     auto element = this->_element->extract_sub_element(component);
@@ -86,17 +91,13 @@ class FunctionSpace
         = std::make_shared<DofMap>(_dofmap->extract_sub_dofmap(component));
 
     // Create new sub space
-    auto sub_space = std::make_shared<FunctionSpace<T>>(_mesh, element, dofmap);
+    FunctionSpace sub_space(_mesh, element, dofmap);
 
     // Set root space id and component w.r.t. root
-    sub_space->_root_space_id = _root_space_id;
-    sub_space->_component = _component;
-    sub_space->_component.insert(sub_space->_component.end(), component.begin(),
-                                 component.end());
-
-    // Insert new subspace into cache
-    _subspaces.emplace(sub_space->_component, sub_space);
-
+    sub_space._root_space_id = _root_space_id;
+    sub_space._component = _component;
+    sub_space._component.insert(sub_space._component.end(), component.begin(),
+                                component.end());
     return sub_space;
   }
 
@@ -138,7 +139,7 @@ class FunctionSpace
   /// Collapse a subspace and return a new function space and a map from
   /// new to old dofs
   /// @return The new function space and a map from new to old dofs
-  std::pair<FunctionSpace<T>, std::vector<std::int32_t>> collapse() const
+  std::pair<FunctionSpace, std::vector<std::int32_t>> collapse() const
   {
     if (_component.empty())
       throw std::runtime_error("Function space is not a subspace");
@@ -169,7 +170,7 @@ class FunctionSpace
   /// @return The dof coordinates `[([x0, y0, z0], [x1, y1, z1], ...)`
   /// if `transpose` is false, and otherwise the returned data is
   /// transposed. Storage is row-major.
-  std::vector<T> tabulate_dof_coordinates(bool transpose) const
+  std::vector<geometry_type> tabulate_dof_coordinates(bool transpose) const
   {
     if (!_component.empty())
     {
@@ -215,28 +216,30 @@ class FunctionSpace
     // Get coordinate map
     if (_mesh->geometry().cmaps().size() > 1)
       throw std::runtime_error("Mixed topology not supported");
-    const CoordinateElement<T>& cmap = _mesh->geometry().cmaps()[0];
+    const CoordinateElement<geometry_type>& cmap = _mesh->geometry().cmaps()[0];
 
     // Prepare cell geometry
     auto x_dofmap = _mesh->geometry().dofmap();
     const std::size_t num_dofs_g = cmap.dim();
-    std::span<const T> x_g = _mesh->geometry().x();
+    std::span<const geometry_type> x_g = _mesh->geometry().x();
 
     // Array to hold coordinates to return
     const std::size_t shape_c0 = transpose ? 3 : num_dofs;
     const std::size_t shape_c1 = transpose ? num_dofs : 3;
-    std::vector<T> coords(shape_c0 * shape_c1, 0);
+    std::vector<geometry_type> coords(shape_c0 * shape_c1, 0);
 
     using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
-        T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
+        geometry_type,
+        MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
     using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
-        const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
+        const geometry_type,
+        MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
 
     // Loop over cells and tabulate dofs
-    std::vector<T> x_b(scalar_dofs * gdim);
+    std::vector<geometry_type> x_b(scalar_dofs * gdim);
     mdspan2_t x(x_b.data(), scalar_dofs, gdim);
 
-    std::vector<T> coordinate_dofs_b(num_dofs_g * gdim);
+    std::vector<geometry_type> coordinate_dofs_b(num_dofs_g * gdim);
     mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
 
     auto map = _mesh->topology()->index_map(tdim);
@@ -251,11 +254,12 @@ class FunctionSpace
     }
 
     auto apply_dof_transformation
-        = _element->template get_pre_dof_transformation_function<T>();
+        = _element
+              ->template get_pre_dof_transformation_function<geometry_type>();
 
     const std::array<std::size_t, 4> phi_shape
         = cmap.tabulate_shape(0, Xshape[0]);
-    std::vector<T> phi_b(
+    std::vector<geometry_type> phi_b(
         std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
     cmdspan4_t phi_full(phi_b.data(), phi_shape);
     cmap.tabulate(0, X, Xshape, phi_b);
@@ -300,20 +304,26 @@ class FunctionSpace
   }
 
   /// The mesh
-  std::shared_ptr<const mesh::Mesh<T>> mesh() const { return _mesh; }
+  std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const
+  {
+    return _mesh;
+  }
 
   /// The finite element
-  std::shared_ptr<const FiniteElement<T>> element() const { return _element; }
+  std::shared_ptr<const FiniteElement<geometry_type>> element() const
+  {
+    return _element;
+  }
 
   /// The dofmap
   std::shared_ptr<const DofMap> dofmap() const { return _dofmap; }
 
 private:
   // The mesh
-  std::shared_ptr<const mesh::Mesh<T>> _mesh;
+  std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
 
   // The finite element
-  std::shared_ptr<const FiniteElement<T>> _element;
+  std::shared_ptr<const FiniteElement<geometry_type>> _element;
 
   // The dofmap
   std::shared_ptr<const DofMap> _dofmap;
@@ -324,9 +334,6 @@ class FunctionSpace
   // Unique identifier for the space and for its root space
   boost::uuids::uuid _id;
   boost::uuids::uuid _root_space_id;
-
-  // Cache of subspaces
-  mutable std::map<std::vector<int>, std::weak_ptr<FunctionSpace>> _subspaces;
 };
 
 /// Extract FunctionSpaces for (0) rows blocks and (1) columns blocks