From 30d09ec108afaf7e02d0cac207c805dbaed0aad4 Mon Sep 17 00:00:00 2001 From: Igor Baratta Date: Fri, 1 Dec 2023 13:07:42 +0000 Subject: [PATCH] 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 --- cpp/demo/hyperelasticity/main.cpp | 2 +- cpp/dolfinx/fem/Constant.h | 11 ++-- cpp/dolfinx/fem/FiniteElement.h | 31 ++++++----- cpp/dolfinx/fem/Function.h | 5 +- cpp/dolfinx/fem/FunctionSpace.h | 89 +++++++++++++++++-------------- python/demo/demo_mixed-poisson.py | 12 +++-- python/demo/demo_stokes.py | 17 +++--- python/dolfinx/fem/function.py | 24 ++++++--- 8 files changed, 110 insertions(+), 81 deletions(-) diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp index 8517c3e6a4d..f4f675a4a9c 100644 --- 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>(std::vector{0, 0, 0}, bdofs_left, V), std::make_shared>(u_rotation, bdofs_right)}; diff --git a/cpp/dolfinx/fem/Constant.h b/cpp/dolfinx/fem/Constant.h index 2095c68d9d7..a344f1787e1 100644 --- a/cpp/dolfinx/fem/Constant.h +++ b/cpp/dolfinx/fem/Constant.h @@ -22,13 +22,16 @@ template 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 c) + explicit Constant(std::span c) : Constant(c, std::vector{c.size()}) { } @@ -36,13 +39,13 @@ class Constant /// @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 c, std::span shape) + Constant(std::span c, std::span shape) : value(c.begin(), c.end()), shape(shape.begin(), shape.end()) { } /// Values, stored as a row-major flattened array - std::vector value; + std::vector value; /// Shape std::vector shape; diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h index c69c904c1b0..9433b605f7b 100644 --- a/cpp/dolfinx/fem/FiniteElement.h +++ b/cpp/dolfinx/fem/FiniteElement.h @@ -26,6 +26,9 @@ template 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& element, + FiniteElement(const basix::FiniteElement& element, const std::vector& 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 values, std::span X, + void tabulate(std::span values, + std::span X, std::array 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::array> - tabulate(std::span X, std::array shape, + std::pair, std::array> + tabulate(std::span X, std::array 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>>& + const std::vector>>& sub_elements() const noexcept; /// Extract sub finite element for component - std::shared_ptr> + std::shared_ptr> extract_sub_element(const std::vector& component) const; /// Return underlying basix element (if it exists) - const basix::FiniteElement& basix_element() const; + const basix::FiniteElement& 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::array> + std::pair, std::array> 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::array> + std::pair, std::array> 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::array> + std::pair, std::array> 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>> _sub_elements; + std::vector>> + _sub_elements; // Dimension of each value space std::vector _value_shape; @@ -685,10 +690,10 @@ class FiniteElement bool _needs_dof_transformations; // Basix Element (nullptr for mixed elements) - std::unique_ptr> _element; + std::unique_ptr> _element; // Quadrature points of a quadrature element (0 dimensional array for // all elements except quadrature elements) - std::pair, std::array> _points; + std::pair, std::array> _points; }; } // namespace dolfinx::fem diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h index c8758d3fed4..8c246454b63 100644 --- 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>( + _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 index 76d710f61d5..1ba6cb9340b 100644 --- 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 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> mesh, - std::shared_ptr> element, + FunctionSpace(std::shared_ptr> mesh, + std::shared_ptr> element, std::shared_ptr 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> sub(const std::vector& 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& 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->extract_sub_dofmap(component)); // Create new sub space - auto sub_space = std::make_shared>(_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, std::vector> collapse() const + std::pair> 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 tabulate_dof_coordinates(bool transpose) const + std::vector 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& cmap = _mesh->geometry().cmaps()[0]; + const CoordinateElement& 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 x_g = _mesh->geometry().x(); + std::span 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 coords(shape_c0 * shape_c1, 0); + std::vector coords(shape_c0 * shape_c1, 0); using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + geometry_type, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< - const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + const geometry_type, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; // Loop over cells and tabulate dofs - std::vector x_b(scalar_dofs * gdim); + std::vector x_b(scalar_dofs * gdim); mdspan2_t x(x_b.data(), scalar_dofs, gdim); - std::vector coordinate_dofs_b(num_dofs_g * gdim); + std::vector 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(); + = _element + ->template get_pre_dof_transformation_function(); const std::array phi_shape = cmap.tabulate_shape(0, Xshape[0]); - std::vector phi_b( + std::vector 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> mesh() const { return _mesh; } + std::shared_ptr> mesh() const + { + return _mesh; + } /// The finite element - std::shared_ptr> element() const { return _element; } + std::shared_ptr> element() const + { + return _element; + } /// The dofmap std::shared_ptr dofmap() const { return _dofmap; } private: // The mesh - std::shared_ptr> _mesh; + std::shared_ptr> _mesh; // The finite element - std::shared_ptr> _element; + std::shared_ptr> _element; // The dofmap std::shared_ptr _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::weak_ptr> _subspaces; }; /// Extract FunctionSpaces for (0) rows blocks and (1) columns blocks diff --git a/python/demo/demo_mixed-poisson.py b/python/demo/demo_mixed-poisson.py index b5434288c66..5458e3561ed 100644 --- a/python/demo/demo_mixed-poisson.py +++ b/python/demo/demo_mixed-poisson.py @@ -113,11 +113,13 @@ a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx L = -inner(f, v) * dx +# Get subspace of V +V0 = V.sub(0) fdim = domain.topology.dim - 1 facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0)) -Q, _ = V.sub(0).collapse() -dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top) +Q, _ = V0.collapse() +dofs_top = fem.locate_dofs_topological((V0, Q), fdim, facets_top) def f1(x): @@ -128,11 +130,11 @@ def f1(x): f_h1 = fem.Function(Q) f_h1.interpolate(f1) -bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0)) +bc_top = fem.dirichletbc(f_h1, dofs_top, V0) facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0)) -dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom) +dofs_bottom = fem.locate_dofs_topological((V0, Q), fdim, facets_bottom) def f2(x): @@ -143,7 +145,7 @@ def f2(x): f_h2 = fem.Function(Q) f_h2.interpolate(f2) -bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0)) +bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V0) bcs = [bc_top, bc_bottom] diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py index dcfca2e8916..6cb21fee7f1 100644 --- a/python/demo/demo_stokes.py +++ b/python/demo/demo_stokes.py @@ -460,18 +460,19 @@ def mixed_direct(): W = functionspace(msh, TH) # No slip boundary condition - W0, _ = W.sub(0).collapse() - noslip = Function(W0) + W0 = W.sub(0) + Q, _ = W0.collapse() + noslip = Function(Q) facets = locate_entities_boundary(msh, 1, noslip_boundary) - dofs = locate_dofs_topological((W.sub(0), W0), 1, facets) - bc0 = dirichletbc(noslip, dofs, W.sub(0)) + dofs = locate_dofs_topological((W0, Q), 1, facets) + bc0 = dirichletbc(noslip, dofs, W0) # Driving velocity condition u = (1, 0) on top boundary (y = 1) - lid_velocity = Function(W0) + lid_velocity = Function(Q) lid_velocity.interpolate(lid_velocity_expression) facets = locate_entities_boundary(msh, 1, lid) - dofs = locate_dofs_topological((W.sub(0), W0), 1, facets) - bc1 = dirichletbc(lid_velocity, dofs, W.sub(0)) + dofs = locate_dofs_topological((W0, Q), 1, facets) + bc1 = dirichletbc(lid_velocity, dofs, W0) # Collect Dirichlet boundary conditions bcs = [bc0, bc1] @@ -479,7 +480,7 @@ def mixed_direct(): # Define variational problem (u, p) = ufl.TrialFunctions(W) (v, q) = ufl.TestFunctions(W) - f = Function(W0) + f = Function(Q) a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx) L = form(inner(f, v) * dx) diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 5efb2923683..827afc6aa96 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -447,15 +447,21 @@ def __str__(self): return self.name def sub(self, i: int) -> Function: - """Return a sub-function. + """Return a sub-function (a view into the `Function`). + + Sub-functions are indexed `i = 0, ..., N-1`, where `N` is the + number of sub-spaces. Args: i: Index of the sub-function to extract. - Note: - The sub-functions are numbered i = 0, ..., N-1, where N is - the number of sub-spaces. + Returns: + A view into the parent `Function`. + Note: + If the sub-Function is re-used, for performance reasons the + returned `Function` should be stored by the caller to avoid + repeated re-computation of the subspac. """ return Function(self._V.sub(i), self.x, name=f"{str(self)}_{i}") @@ -566,7 +572,7 @@ def __init__(self, mesh: Mesh, element: ufl.FiniteElementBase, Note: This initialiser is for internal use and not normally called - in user code. Use :func:`FunctionSpace` to create a function space. + in user code. Use :func:`functionspace` to create a function space. Args: mesh: Mesh that space is defined on @@ -615,11 +621,15 @@ def sub(self, i: int) -> FunctionSpace: """Return the i-th sub space. Args: - i: The subspace index + i: Index of the subspace to extract. Returns: - A subspace + A subspace. + 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. """ assert self.ufl_element().num_sub_elements > i sub_element = self.ufl_element().sub_elements[i]