Document shapes of xtensors

cpp/basix/cell.h

@@ -29,7 +29,7 @@ enum class type
 
 /// Cell geometry
 /// @param celltype Cell Type
-/// @return Set of vertex points of the cell
+/// @return Set of vertex points of the cell. Shape is (npoints, gdim)
 xt::xtensor<double, 2> geometry(cell::type celltype);
 
 /// Cell topology
@@ -54,7 +54,7 @@ sub_entity_connectivity(cell::type celltype);
 /// @param celltype The cell::type
 /// @param dim Dimension of sub-entity
 /// @param index Local index of sub-entity
-/// @return Set of vertex points of the sub-entity
+/// @return Set of vertex points of the sub-entity. Shape is (npoints, gdim)
 xt::xtensor<double, 2> sub_entity_geometry(cell::type celltype, int dim,
                                            int index);
 
@@ -86,13 +86,13 @@ double volume(cell::type cell_type);
 
 /// Get the (outward) normals to the facets of a reference cell
 /// @param cell_type Type of cell
-/// @return The outward normals
+/// @return The outward normals. Shape is (nfacets, gdim)
 xt::xtensor<double, 2> facet_outward_normals(cell::type cell_type);
 
 /// Get the normals to the facets of a reference cell oriented using the
 /// low-to-high ordering of the facet
 /// @param cell_type Type of cell
-/// @return The normals
+/// @return The normals. Shape is (nfacets, gdim)
 xt::xtensor<double, 2> facet_normals(cell::type cell_type);
 
 /// Get an array of bools indicating whether or not the facet normals are
@@ -104,16 +104,17 @@ std::vector<bool> facet_orientations(cell::type cell_type);
 /// Get the reference volumes of the facets of a reference cell
 /// @param cell_type Type of cell
 /// @return The volumes of the references associated with each facet
+// FIXME: This could be a std::vector<double>
 xt::xtensor<double, 1> facet_reference_volumes(cell::type cell_type);
 
 /// Get the types of the subentities of a reference cell
 /// @param cell_type Type of cell
-/// @return The subentity types
+/// @return The subentity types. Indices are (tdim, entity)
 std::vector<std::vector<cell::type>> subentity_types(cell::type cell_type);
 
 /// Get the jacobians of the facets of a reference cell
 /// @param cell_type Type of cell
-/// @return The jacobians of the facets
+/// @return The jacobians of the facets. Shape is (nfacets, gdim, gdim - 1)
 xt::xtensor<double, 3> facet_jacobians(cell::type cell_type);
 
 } // namespace basix::cell

cpp/basix/dof-transformations.h

@@ -17,13 +17,18 @@ namespace basix::doftransforms
 
 /// Compute the entity DOF transformations for an element
 /// @param[in] cell_type The cell type
-/// @param[in] x Interpolation points for the element
-/// @param[in] M Interpolation matrix fot the element
+/// @param[in] x Interpolation points for the element. Indices are (tdim, entity
+/// index, point index, dim)
+/// @param[in] M Interpolation matrix for the element. Indices are (tdim, entity
+/// index, dof, vs, point_index, derivative)
 /// @param[in] coeffs The coefficients that define the basis functions of the
-/// element in terms of the orthonormal basis
+/// element in terms of the orthonormal basis. Shape is (dim(Legendre
+/// polynomials), dim(finite element polyset))
 /// @param[in] degree The degree of the element
 /// @param[in] vs The value size of the element
 /// @param[in] map_type The map type used by the element
+/// @return Entity transformations. For each cell, the shape is
+/// (ntransformation, ndofs, ndofs)
 std::map<cell::type, xt::xtensor<double, 3>> compute_entity_transformations(
     cell::type cell_type,
     const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,

cpp/basix/finite-element.h

@@ -29,13 +29,15 @@ namespace element
 /// version of the element. This discontinuous version will have the
 /// same DOFs but they will all be associated with the interior of the
 /// reference cell.
-/// @param[in] x Interpolation points. Shape is (tdim, entity index,
+/// @param[in] x Interpolation points. Indices are (tdim, entity index,
 /// point index, dim)
 /// @param[in] M The interpolation matrices. Indices are (tdim, entity
 /// index, dof, vs, point_index, derivative)
 /// @param[in] tdim The topological dimension of the cell the element is
 /// defined on
 /// @param[in] value_size The value size of the element
+/// @return Versions of x and M that define a discontinuous version of the
+/// element (with the same shapes as x and M)
 std::tuple<std::array<std::vector<xt::xtensor<double, 2>>, 4>,
            std::array<std::vector<xt::xtensor<double, 4>>, 4>>
 make_discontinuous(const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
@@ -202,8 +204,9 @@ class FiniteElement
   /// @param[in] value_shape The value shape of the element
   /// @param[in] wcoeffs Matrices for the kth value index containing the
   /// expansion coefficients defining a polynomial basis spanning the
-  /// polynomial space for this element
-  /// @param[in] x Interpolation points. Shape is (tdim, entity index,
+  /// polynomial space for this element. Shape is (dim(Legendre polynomials),
+  /// dim(finite element polyset))
+  /// @param[in] x Interpolation points. Indices are (tdim, entity index,
   /// point index, dim)
   /// @param[in] M The interpolation matrices. Indices are (tdim, entity
   /// index, dof, vs, point_index, derivative)
@@ -440,7 +443,8 @@ class FiniteElement
   /// @param[in] J The Jacobian of the mapping
   /// @param[in] detJ The determinant of the Jacobian of the mapping
   /// @param[in] K The inverse of the Jacobian of the mapping
-  /// @return The function values on the reference
+  /// @return The function values on the reference. The indices are [Jacobian
+  /// index, point index, components].
   xt::xtensor<double, 3> pull_back(const xt::xtensor<double, 3>& u,
                                    const xt::xtensor<double, 3>& J,
                                    const xtl::span<const double>& detJ,
@@ -619,11 +623,13 @@ class FiniteElement
   ///   reflection: [[0, 1],
   ///                [1, 0]]
   /// ~~~~~~~~~~~~~~~~
-  /// @return The base transformations for this element
+  /// @return The base transformations for this element. The shape is
+  /// (ntranformations, ndofs, ndofs)
   xt::xtensor<double, 3> base_transformations() const;
 
   /// Return the entity dof transformation matrices
-  /// @return The entity transformations for the subentities of this element
+  /// @return The entity transformations for the subentities of this element.
+  /// The shape for each cell is (ntransformations, ndofs, ndofs)
   std::map<cell::type, xt::xtensor<double, 3>> entity_transformations() const;
 
   /// Permute the dof numbering on a cell
@@ -799,14 +805,15 @@ class FiniteElement
   ///     coefficients[::block_size] = i_m * values
   /// \endcode
   ///
-  /// @return The interpolation matrix
+  /// @return The interpolation matrix. Shape is (ndofs, number of interpolation
+  /// points)
   const xt::xtensor<double, 2>& interpolation_matrix() const;
 
   /// Get the dual matrix.
   ///
   /// This is the matrix @f$BD^{T}@f$, as described in the documentation
   /// of the `FiniteElement()` constructor.
-  /// @return The dual matrix
+  /// @return The dual matrix. Shape is (ndofs, ndofs)
   const xt::xtensor<double, 2>& dual_matrix() const;
 
   /// Get the coefficients that define the polynomial set in terms of the
@@ -840,11 +847,13 @@ class FiniteElement
   ///
   /// These coefficients are only stored for custom elements. This function will
   /// throw an exception if called on a non-custom element
+  /// @return Coefficient matrix. Shape is (dim(Lagrange polynomials),
+  /// dim(finite element polyset))
   const xt::xtensor<double, 2>& wcoeffs() const;
 
   /// Get the interpolation points for each subentity.
   ///
-  /// The shape of this data is (tdim, entity index, point index, dim).
+  /// The indices of this data are (tdim, entity index, point index, dim).
   const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x() const;
 
   /// Get the interpolation matrices for each subentity.
@@ -881,13 +890,15 @@ class FiniteElement
   ///
   /// These matrices are only stored for custom elements. This function will
   /// throw an exception if called on a non-custom element
+  /// @return The interpolation matrices. The indices of this data are (tdim,
+  /// entity index, dof, vs, point_index, derivative)
   const std::array<std::vector<xt::xtensor<double, 4>>, 4>& M() const;
 
   /// Get the matrix of coefficients.
   ///
   /// This is the matrix @f$C@f$, as described in the documentation of
   /// the `FiniteElement()` constructor.
-  /// @return The dual matrix
+  /// @return The coefficient matrix. Shape is (ndofs, ndofs)
   const xt::xtensor<double, 2>& coefficient_matrix() const;
 
 
@@ -1067,11 +1078,12 @@ class FiniteElement
 /// @param[in] value_shape The value shape of the element
 /// @param[in] wcoeffs Matrices for the kth value index containing the
 /// expansion coefficients defining a polynomial basis spanning the
-/// polynomial space for this element
-/// @param[in] x Interpolation points. Shape is (tdim, entity index,
+/// polynomial space for this element. Shape is (dim(Legendre polynomials),
+/// dim(finite element polyset))
+/// @param[in] x Interpolation points. Indices are (tdim, entity index,
 /// point index, dim)
 /// @param[in] M The interpolation matrices. Indices are (tdim, entity
-/// index, dof, vs, point_index)
+/// index, dof, vs, point_index, derivative)
 /// @param[in] interpolation_nderivs The number of derivatives that need to be
 /// used during interpolation
 /// @param[in] map_type The type of map to be used to map values from

cpp/basix/interpolation.h

@@ -38,7 +38,8 @@ class FiniteElement;
 /// @param[in] element_from The element to interpolate from
 /// @param[in] element_to The element to interpolate to
 /// @return Matrix operator that maps the 'from' degrees-of-freedom to
-/// the 'to' degrees-of-freedom
+/// the 'to' degrees-of-freedom. Shape is (ndofs(element_to),
+/// ndofs(element_from))
 xt::xtensor<double, 2>
 compute_interpolation_operator(const FiniteElement& element_from,
                                const FiniteElement& element_to);

cpp/basix/lattice.h

@@ -82,7 +82,7 @@ enum class simplex_method
 /// @param type A lattice type
 /// @param exterior If set, includes outer boundaries
 /// @param simplex_method The method used to generate points on simplices
-/// @return Set of points
+/// @return Set of points. Shape is (npoints, tdim)
 xt::xtensor<double, 2>
 create(cell::type celltype, int n, lattice::type type, bool exterior,
        lattice::simplex_method simplex_method = lattice::simplex_method::none);

cpp/basix/moments.h

@@ -31,7 +31,10 @@ namespace moments
 /// being defined
 /// @param value_size The value size of the space being defined
 /// @param q_deg The quadrature degree used for the integrals
-/// @return (interpolation points, interpolation matrix)
+/// @return (interpolation points, interpolation matrix). The indices of the
+/// interpolation points are (number of entities, npoints, gdim). The indices on
+/// the interpolation matrix are (number of entities, ndofs, value_size,
+/// npoints, derivative)
 std::pair<std::vector<xt::xtensor<double, 2>>,
           std::vector<xt::xtensor<double, 4>>>
 make_integral_moments(const FiniteElement& moment_space, cell::type celltype,
@@ -54,7 +57,10 @@ make_integral_moments(const FiniteElement& moment_space, cell::type celltype,
 /// defined
 /// @param value_size The value size of the space being defined
 /// @param q_deg The quadrature degree used for the integrals
-/// @return (interpolation points, interpolation matrix)
+/// @return (interpolation points, interpolation matrix). The indices of the
+/// interpolation points are (number of entities, npoints, gdim). The indices on
+/// the interpolation matrix are (number of entities, ndofs, value_size,
+/// npoints, derivative)
 std::pair<std::vector<xt::xtensor<double, 2>>,
           std::vector<xt::xtensor<double, 4>>>
 make_dot_integral_moments(const FiniteElement& V, cell::type celltype,
@@ -71,7 +77,10 @@ make_dot_integral_moments(const FiniteElement& V, cell::type celltype,
 /// @param value_size The value size of the space being defined the
 /// space
 /// @param q_deg The quadrature degree used for the integrals
-/// @return (interpolation points, interpolation matrix)
+/// @return (interpolation points, interpolation matrix). The indices of the
+/// interpolation points are (number of entities, npoints, gdim). The indices on
+/// the interpolation matrix are (number of entities, ndofs, value_size,
+/// npoints, derivative)
 std::pair<std::vector<xt::xtensor<double, 2>>,
           std::vector<xt::xtensor<double, 4>>>
 make_tangent_integral_moments(const FiniteElement& V, cell::type celltype,
@@ -87,7 +96,10 @@ make_tangent_integral_moments(const FiniteElement& V, cell::type celltype,
 /// being defined
 /// @param[in] value_size The value size of the space being defined
 /// @param[in] q_deg The quadrature degree used for the integrals
-/// @return (interpolation points, interpolation matrix)
+/// @return (interpolation points, interpolation matrix). The indices of the
+/// interpolation points are (number of entities, npoints, gdim). The indices on
+/// the interpolation matrix are (number of entities, ndofs, value_size,
+/// npoints, derivative)
 std::pair<std::vector<xt::xtensor<double, 2>>,
           std::vector<xt::xtensor<double, 4>>>
 make_normal_integral_moments(const FiniteElement& V, cell::type celltype,

cpp/basix/quadrature.h

@@ -54,12 +54,12 @@ quadrature::type get_default_rule(cell::type celltype, int m);
 
 /// Get Gauss-Lobatto-Legendre (GLL) points on the interval [0, 1].
 /// @param[in] m The number of points
-/// @return An array of GLL points
+/// @return An array of GLL points. Shape is (num points, gdim)
 xt::xtensor<double, 2> get_gll_points(int m);
 
 /// Get Gauss-Legendre (GL) points on the interval [0, 1].
 /// @param[in] m The number of points
-/// @return An array of GL points
+/// @return An array of GL points. Shape is (num points, gdim)
 xt::xtensor<double, 2> get_gl_points(int m);
 
 } // namespace basix::quadrature

python/docs.h

@@ -20,7 +20,7 @@ Cell geometry
     celltype (basix.CellType): Cell Type
 
 Returns::
-    numpy.ndarray[numpy.float64]: Set of vertex points of the cell
+    numpy.ndarray[numpy.float64]: Set of vertex points of the cell. Shape is (npoints, gdim)
 )";
 
 const std::string sub_entity_connectivity = R"(
@@ -48,7 +48,7 @@ Sub-entity of a cell, given by topological dimension and index
     index (int): Local index of sub-entity
 
 Returns:
-    List[List[List[List[int]]]]: Set of vertex points of the sub-entity
+    List[List[List[List[int]]]]: Set of vertex points of the sub-entity. Shape is (npoints, gdim)
 )";
 
 const std::string create_lattice__celltype_n_type_exterior = R"(
@@ -72,7 +72,7 @@ the Gauss-Lobatto-Legendre quadrature points. These are the same as
     exterior (bool): If set, includes outer boundaries
 
 Returns:
-    numpy.ndarray[numpy.float64]: Set of points
+    numpy.ndarray[numpy.float64]: Set of points. Shape is (npoints, tdim)
 )";
 
 const std::string create_lattice__celltype_n_type_exterior_method = R"(
@@ -97,7 +97,7 @@ the Gauss-Lobatto-Legendre quadrature points. These are the same as
     simplex_method (basix.LatticeSimplexMethod): The method used to generate points on simplices
 
 Returns:
-    numpy.ndarray[numpy.float64]: Set of points
+    numpy.ndarray[numpy.float64]: Set of points. Shape is (npoints, tdim)
 )";
 
 const std::string cell_volume = R"(
@@ -118,7 +118,7 @@ low-to-high ordering of the facet
     cell_type (basix.CellType): Type of cell
 
 Returns:
-    numpy.ndarray[numpy.float64]: The normals
+    numpy.ndarray[numpy.float64]: The normals. Shape is (nfacets, gdim)
 )";
 
 const std::string cell_facet_reference_volumes = R"(
@@ -138,7 +138,7 @@ Get the (outward) normals to the facets of a reference cell
     cell_type (basix.CellType): Type of cell
 
 Returns:
-    numpy.ndarray[numpy.float64]: The outward normals
+    numpy.ndarray[numpy.float64]: The outward normals. Shape is (nfacets, gdim)
 )";
 
 const std::string cell_facet_orientations = R"(
@@ -159,7 +159,7 @@ Get the jacobians of the facets of a reference cell
     cell_type (basix.CellType): Type of cell
 
 Returns:
-    numpy.ndarray[numpy.float64]: The jacobians of the facets
+    numpy.ndarray[numpy.float64]: The jacobians of the facets. Shape is (nfacets, gdim, gdim - 1)
 )";
 
 const std::string FiniteElement__tabulate = R"(
@@ -213,7 +213,8 @@ Map function values from a physical cell to the reference
     K (numpy.ndarray[numpy.float64]): The inverse of the Jacobian of the mapping
 
 Returns:
-    numpy.ndarray[numpy.float64]: The function values on the reference
+    numpy.ndarray[numpy.float64]: The function values on the reference. The indices are [Jacobian
+    index, point index, components].
 )";
 
 const std::string FiniteElement__apply_dof_transformation = R"(
@@ -358,14 +359,16 @@ For these DOFs, the subblocks of the base transformation matrices are:
 
 
 Returns:
-    numpy.ndarray[numpy.float64]: The base transformations for this element
+    numpy.ndarray[numpy.float64]: The base transformations for this element. The shape is
+    (ntranformations, ndofs, ndofs)
 )";
 
 const std::string FiniteElement__entity_transformations = R"(
 Return the entity dof transformation matrices
 
 Returns:
-    dict: The base transformations for this element
+    dict: The base transformations for this element. The shape is
+    (ntranformations, ndofs, ndofs)
 )";
 
 const std::string FiniteElement__get_tensor_product_representation = R"(
@@ -390,9 +393,9 @@ Create a custom finite element
 Args:
     cell_type (basix.CellType): The cell type
     value_shape (List[int]): The value shape of the element
-    wcoeffs (numpy.ndarray[numpy.float64]): Matrices for the kth value index containing the expansion coefficients defining a polynomial basis spanning the polynomial space for this element
-    x (List[List[numpy.ndarray[numpy.float64]]]): Interpolation points. Shape is (tdim, entity index, point index, dim)
-    M (List[List[numpy.ndarray[numpy.float64]]]): The interpolation matrices. Indices are (tdim, entity index, dof, vs, point_index)
+    wcoeffs (numpy.ndarray[numpy.float64]): Matrices for the kth value index containing the expansion coefficients defining a polynomial basis spanning the polynomial space for this element. Shape is (dim(Legendre polynomials), dim(finite element polyset))
+    x (List[List[numpy.ndarray[numpy.float64]]]): Interpolation points. Indices are (tdim, entity index, point index, dim)
+    M (List[List[numpy.ndarray[numpy.float64]]]): The interpolation matrices. Indices are (tdim, entity index, dof, vs, point_index, derivative)
     interpolation_nderivs (int): The number of derivatives that need to be used during interpolation
     map_type (basix.MapType): The type of map to be used to map values from the reference to a cell
     discontinuous (bool): Indicates whether or not this is the discontinuous version of the element
@@ -548,7 +551,8 @@ each other out.
 
 Returns:
     numpy.ndarray[numpy.float64]: Matrix operator that maps the 'from' degrees-of-freedom to
-    the 'to' degrees-of-freedom
+    the 'to' degrees-of-freedom. Shape is (ndofs(element_to),
+    ndofs(element_from))
 )";
 
 const std::string tabulate_polynomial_set = R"(