Add iso macroelement on an interval (#680)

* Add polyset::type input to all polyset functions, pass standard in everywhere for now

* Add polyset_type input to FiniteElement constructor

* implement quadrature for macro element on interval

* implement macro polyset on an interval

* Add test, update demos, clang format

* update _basixcpp.pyi

* tabulate_polynomial_set

* Remove overload of make_quadrature, add python function with optional args instead

* default value

* flake8

* correct _basixcpp.pyi

* rule=

* update demos

* flake8

* flake8

* Add new element to test_create

* update README

* dolfinx branch

* Add polyset_type input to custom element

* flake

* _

* implement polyset_type for more elements

* update tests

* pt -> ptype

* update _basixcpp.pyi

* update demos

* document arg

* add "iso" to family strings

* _basixcpp.pyi

* move basix.quadrature.make_quadrature back to basix.make_quadrature

Author: mscroggs

.github/workflows/dolfin-tests.yml

@@ -51,7 +51,7 @@ jobs:
         if: github.event_name != 'workflow_dispatch'
         run: |
           python3 -m pip install git+https://github.com/FEniCS/ufl.git
-          python3 -m pip install git+https://github.com/FEniCS/ffcx.git
+          python3 -m pip install git+https://github.com/FEniCS/ffcx.git@mscroggs/make_quadrature
       - name: Install FEniCS Python components
         if: github.event_name == 'workflow_dispatch'
         run: |
@@ -64,7 +64,7 @@ jobs:
         with:
           path: ./dolfinx
           repository: FEniCS/dolfinx
-          ref: main
+          ref: mscroggs/make_quadrature
       - name: Get DOLFINx
         if: github.event_name == 'workflow_dispatch'
         uses: actions/checkout@v3

.github/workflows/ffcx-tests.yml

@@ -54,7 +54,7 @@ jobs:
         with:
           path: ./ffcx
           repository: FEniCS/ffcx
-          ref: main
+          ref: mscroggs/make_quadrature
       - name: Get FFCx source (specified branch)
         if: github.event_name == 'workflow_dispatch'
         uses: actions/checkout@v3

README.md

@@ -71,7 +71,7 @@ The following elements are supported on an interval:
   - [Bubble](https://defelement.com/elements/bubble.html)
   - [Serendipity](https://defelement.com/elements/serendipity.html)
   - [Hermite](https://defelement.com/elements/hermite.html)
-
+  - [iso](https://defelement.com/elements/p1-iso-p2.html)
 
 ### Triangle
 

cpp/basix/docs.h

@@ -54,7 +54,7 @@ Sub-entity of a cell, given by topological dimension and index
 )";
 
 const std::string create_lattice__celltype_n_type_exterior = R"(
-@brief Create a lattice of points on a reference cell optionally
+Create a lattice of points on a reference cell optionally
 including the outer surface points.
 
 For a given `celltype`, this creates a set of points on a regular
@@ -80,7 +80,7 @@ type::equispaced when `n < 3`.
 )";
 
 const std::string create_lattice__celltype_n_type_exterior_method = R"(
-@brief Create a lattice of points on a reference cell optionally
+Create a lattice of points on a reference cell optionally
 including the outer surface points.
 
 For a given `celltype`, this creates a set of points on a regular
@@ -169,7 +169,7 @@ Get the jacobians of the facets of a reference cell
 )";
 
 const std::string FiniteElement__tabulate = R"(
-@brief Compute basis values and derivatives at set of points.
+Compute basis values and derivatives at set of points.
 
 NOTE: The version of `FiniteElement::tabulate` with the basis data
 as an out argument should be preferred for repeated call where
@@ -270,7 +270,7 @@ performance is critical.
 )";
 
 const std::string FiniteElement__base_transformations = R"(
-@brief Get the base transformations.
+Get the base transformations.
 
 The base transformations represent the effect of rotating or reflecting
 a subentity of the cell on the numbering and orientation of the DOFs.
@@ -409,6 +409,7 @@ Create a custom finite element
     discontinuous: Indicates whether or not this is the discontinuous version of the element
     highest_complete_degree: The highest degree n such that a Lagrange (or vector Lagrange) element of degree n is a subspace of this element
     highest_degree: The degree of a polynomial in this element's polyset
+    poly_type: The type of polyset to use for this element
 
 Returns:
     A custom finite element
@@ -434,30 +435,30 @@ Create an element using a given Lagrange variant and a given DPC variant
 )";
 
 const std::string compute_interpolation_operator = R"(
-Computes a matrix that represents the interpolation between two
-elements.
+Compute a matrix that represents the interpolation between
+two elements.
 
 If the two elements have the same value size, this function returns
 the interpolation between them.
 
-If element_from has value size 1 and element_to has value size > 1, then
-this function returns a matrix to interpolate from a blocked element_from
-(ie multiple copies of element_from) into element_to.
+If element_from has value size 1 and element_to has value size > 1,
+then this function returns a matrix to interpolate from a blocked
+element_from (ie multiple copies of element_from) into element_to.
 
-If element_to has value size 1 and element_from has value size > 1, then
-this function returns a matrix that interpolates the components of
-element_from into copies of element_to.
+If element_to has value size 1 and element_from has value size > 1,
+then this function returns a matrix that interpolates the components
+of element_from into copies of element_to.
 
 NOTE: If the elements have different value sizes and both are
 greater than 1, this function throws a runtime error
 
-In order to interpolate functions between finite element spaces on arbitrary
-cells, the functions must be pulled back to the reference element (this pull
-back includes applying DOF transformations). The matrix that this function
-returns can then be applied, then the result pushed forward to the cell. If
-element_from and element_to have the same map type, then only the DOF
-transformations need to be applied, as the pull back and push forward cancel
-each other out.
+In order to interpolate functions between finite element spaces on
+arbitrary cells, the functions must be pulled back to the reference
+element (this pull back includes applying DOF transformations). The
+matrix that this function returns can then be applied, then the
+result pushed forward to the cell. If element_from and element_to
+have the same map type, then only the DOF transformations need to be
+applied, as the pull back and push forward cancel each other out.
 
 Args:
     element_from: The element to interpolate from
@@ -470,7 +471,7 @@ each other out.
 )";
 
 const std::string tabulate_polynomial_set = R"(
-@brief Tabulate the orthonormal polynomial basis, and derivatives,
+Tabulate the orthonormal polynomial basis, and derivatives,
 at points on the reference cell.
 
 All derivatives up to the given order are computed. If derivatives
@@ -486,6 +487,7 @@ there are `(nderiv + 1)(nderiv + 2)(nderiv + 3)/6`. The ordering is
 
 Args:
     celltype: Cell type
+    ptype: The polynomial type
     d: Polynomial degree
     n: Maximum derivative order. Use n = 0 for the basis only.
     x: Points at which to evaluate the basis. The shape is (number of points, geometric dimension).
@@ -512,7 +514,7 @@ there are `(nderiv + 1)(nderiv + 2)(nderiv + 3)/6`. The ordering is
 )";
 
 const std::string tabulate_polynomials = R"(
-@brief Tabulate a set of polynomials.
+Tabulate a set of polynomials.
 
 Args:
     polytype: Polynomial type
@@ -526,7 +528,7 @@ const std::string tabulate_polynomials = R"(
 )";
 
 const std::string polynomials_dim = R"(
-@brief Dimension of a polynomial space.
+Dimension of a polynomial space.
 
 Args:
     polytype: The polynomial type
@@ -538,24 +540,13 @@ const std::string polynomials_dim = R"(
     polynomial degree @p d
 )";
 
-const std::string make_quadrature__rule_celltype_m = R"(
+const std::string make_quadrature__rule_celltype_polytype_m = R"(
 Make a quadrature rule on a reference cell
 
 Args:
     rule: Type of quadrature rule (or use quadrature::Default)
     celltype: The cell type
-    m: Maximum degree of polynomial that this quadrature rule will integrate exactly
-
-Returns:
-    List of points and list of weights. The number of points
-    arrays has shape (num points, gdim)
-)";
-
-const std::string make_quadrature__celltype_m = R"(
-Make a default quadrature rule on reference cell
-
-Args:
-    celltype: The cell type
+    polytype: The polyset type
     m: Maximum degree of polynomial that this quadrature rule will integrate exactly
 
 Returns:
@@ -574,9 +565,9 @@ Compute trivial indexing in a 1D array (for completeness)
 )";
 
 const std::string index__p_q = R"(
-Compute indexing in a 2D triangular array compressed into a 1D array.
-This can be used to find the index of a derivative returned by
-`FiniteElement::tabulate`. For instance to find d2N/dx2, use
+Compute indexing in a 2D triangular array compressed into a 1D
+array. This can be used to find the index of a derivative returned
+by `FiniteElement::tabulate`. For instance to find d2N/dx2, use
 `FiniteElement::tabulate(2, points)[idx(2, 0)];`
 
 Args:
@@ -588,7 +579,8 @@ This can be used to find the index of a derivative returned by
 )";
 
 const std::string index__p_q_r = R"(
-Compute indexing in a 3D tetrahedral array compressed into a 1D array
+Compute indexing in a 3D tetrahedral array compressed into a 1D
+array
 
 Args:
     p: Index in x
@@ -610,4 +602,32 @@ Get the intersection of two Sobolev spaces.
     Intersection of the spaces
 )";
 
+const std::string superset = R"(
+Get the polyset types that is a superset of two types on the given
+cell
+
+Args:
+    cell: The cell type
+    type1: The first polyset type
+    type2: The second polyset type
+
+Returns::
+    The superset type
+)";
+
+const std::string restriction = R"(
+Get the polyset type that represents the restrictions of a type on a
+subentity
+
+Args:
+    ptype: The polyset type
+    cell: The cell type
+    restriction_cell: The cell type of the subentity
+
+Returns::
+    The restricted polyset type
+)";
+
+
+
 } // namespace basix::docstring

cpp/basix/dof-transformations.cpp

@@ -327,7 +327,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
 
   const std::size_t ndofs = imat.extent(0);
   const std::size_t npts = pts.extent(0);
-  const int psize = polyset::dim(cell_type, degree);
+  const int psize = polyset::dim(cell_type, polyset::type::standard, degree);
 
   std::size_t dofstart = 0;
   for (int d = 0; d < tdim; ++d)
@@ -352,7 +352,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
   }
 
   auto [polyset_vals_b, polyset_shape] = polyset::tabulate(
-      cell_type, degree, 0,
+      cell_type, polyset::type::standard, degree, 0,
       mdspan_t<const T, 2>(mapped_pts.data(), mapped_pts.extents()));
   assert(polyset_shape[0] == 1);
   mdspan_t<const T, 2> polyset_vals(polyset_vals_b.data(), polyset_shape[1],

cpp/basix/e-brezzi-douglas-marini.cpp

@@ -41,7 +41,8 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
       = create_lagrange<T>(facettype, degree, lvariant, true);
   {
     auto [_x, xshape, _M, Mshape] = moments::make_normal_integral_moments<T>(
-        facet_moment_space, celltype, tdim, degree * 2);
+        facet_moment_space, celltype, polyset::type::standard, tdim,
+        degree * 2);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -56,8 +57,8 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
   {
     // Interior integral moment
     auto [_x, xshape, _M, Mshape] = moments::make_dot_integral_moments<T>(
-        create_nedelec<T>(celltype, degree - 1, lvariant, true), celltype, tdim,
-        2 * degree - 1);
+        create_nedelec<T>(celltype, degree - 1, lvariant, true), celltype,
+        polyset::type::standard, tdim, 2 * degree - 1);
     assert(_x.size() == _M.size());
     for (std::size_t i = 0; i < _x.size(); ++i)
     {
@@ -90,11 +91,12 @@ FiniteElement<T> element::create_bdm(cell::type celltype, int degree,
   }
 
   // The number of order (degree) scalar polynomials
-  const std::size_t ndofs = tdim * polyset::dim(celltype, degree);
+  const std::size_t ndofs
+      = tdim * polyset::dim(celltype, polyset::type::standard, degree);
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::HDiv;
   return FiniteElement<T>(
-      family::BDM, celltype, degree, {tdim},
+      family::BDM, celltype, polyset::type::standard, degree, {tdim},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 0, maps::type::contravariantPiola, space, discontinuous, degree,
       degree, lvariant, element::dpc_variant::unset);

cpp/basix/e-bubble.cpp

@@ -67,10 +67,11 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
 
   // Evaluate the expansion polynomials at the quadrature points
   const auto [_pts, wts] = quadrature::make_quadrature<T>(
-      quadrature::type::Default, celltype, 2 * degree);
+      quadrature::type::Default, celltype, polyset::type::standard, 2 * degree);
   impl::mdspan_t<const T, 2> pts(_pts.data(), wts.size(),
                                  _pts.size() / wts.size());
-  const auto [_phi, shape] = polyset::tabulate(celltype, degree, 0, pts);
+  const auto [_phi, shape]
+      = polyset::tabulate(celltype, polyset::type::standard, degree, 0, pts);
   impl::mdspan_t<const T, 3> phi(_phi.data(), shape);
 
   // The number of order (degree) polynomials
@@ -103,7 +104,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   {
   case cell::type::interval:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -115,7 +117,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::triangle:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 3, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 3, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -128,7 +131,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::tetrahedron:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 4, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 4, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -142,7 +146,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::quadrilateral:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -155,7 +160,8 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   }
   case cell::type::hexahedron:
   {
-    const auto [_phi1, shape] = polyset::tabulate(celltype, degree - 2, 0, pts);
+    const auto [_phi1, shape] = polyset::tabulate(
+        celltype, polyset::type::standard, degree - 2, 0, pts);
     impl::mdspan_t<const T, 3> p1(_phi1.data(), shape);
     phi1 = create_phi1(p1, phi1_buffer);
     for (std::size_t i = 0; i < pts.extent(0); ++i)
@@ -185,9 +191,10 @@ FiniteElement<T> basix::element::create_bubble(cell::type celltype, int degree,
   sobolev::space space
       = discontinuous ? sobolev::space::L2 : sobolev::space::H1;
   return FiniteElement<T>(
-      element::family::bubble, celltype, degree, {}, wview, impl::to_mdspan(x),
-      impl::to_mdspan(M), 0, maps::type::identity, space, discontinuous, -1,
-      degree, element::lagrange_variant::unset, element::dpc_variant::unset);
+      element::family::bubble, celltype, polyset::type::standard, degree, {},
+      wview, impl::to_mdspan(x), impl::to_mdspan(M), 0, maps::type::identity,
+      space, discontinuous, -1, degree, element::lagrange_variant::unset,
+      element::dpc_variant::unset);
 }
 //-----------------------------------------------------------------------------
 template FiniteElement<float> element::create_bubble(cell::type, int, bool);

cpp/basix/e-crouzeix-raviart.cpp

@@ -88,7 +88,7 @@ FiniteElement<T> basix::element::create_cr(cell::type celltype, int degree,
   }
 
   return FiniteElement<T>(
-      element::family::CR, celltype, 1, {},
+      element::family::CR, celltype, polyset::type::standard, 1, {},
       impl::mdspan_t<T, 2>(math::eye<T>(ndofs).data(), ndofs, ndofs), xview,
       Mview, 0, maps::type::identity, sobolev::space::L2, discontinuous, degree,
       degree, element::lagrange_variant::unset, element::dpc_variant::unset);