Implement iso element on quads (#683)

Author: mscroggs

README.md

@@ -111,6 +111,7 @@ The following elements are supported on a quadrilateral:
   - [Bubble](https://defelement.com/elements/bubble.html)
   - [DPC](https://defelement.com/elements/dpc.html)
   - [Serendipity](https://defelement.com/elements/serendipity.html)
+  - [iso](https://defelement.com/elements/p1-iso-p2.html)
 
 
 ### Tetrahedron

cpp/basix/dof-transformations.cpp

@@ -317,7 +317,8 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
     mdspan_t<const T, 2> coeffs, const mdarray_t<T, 2>& J, T detJ,
     const mdarray_t<T, 2>& K,
     const std::function<std::array<T, 3>(std::span<const T>)> map_point,
-    int degree, int tdim, int entity, std::size_t vs, const maps::type map_type)
+    int degree, int tdim, int entity, std::size_t vs, const maps::type map_type,
+    const polyset::type ptype)
 {
   if (x[tdim].size() == 0 or x[tdim][entity].extent(0) == 0)
     return {{}, {0, 0}};
@@ -327,7 +328,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, polyset::type::standard, degree);
+  const int psize = polyset::dim(cell_type, ptype, degree);
 
   std::size_t dofstart = 0;
   for (int d = 0; d < tdim; ++d)
@@ -352,7 +353,7 @@ std::pair<std::vector<T>, std::array<std::size_t, 2>> compute_transformation(
   }
 
   auto [polyset_vals_b, polyset_shape] = polyset::tabulate(
-      cell_type, polyset::type::standard, degree, 0,
+      cell_type, ptype, 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],
@@ -424,7 +425,7 @@ doftransforms::compute_entity_transformations(
     std::experimental::mdspan<const T,
                               std::experimental::dextents<std::size_t, 2>>
         coeffs,
-    int degree, std::size_t vs, maps::type map_type)
+    int degree, std::size_t vs, maps::type map_type, polyset::type ptype)
 {
   std::map<cell::type, std::pair<std::vector<T>, std::array<std::size_t, 3>>>
       out;
@@ -441,7 +442,7 @@ doftransforms::compute_entity_transformations(
     {
       auto [t2b, _]
           = compute_transformation(cell_type, x, M, coeffs, J, detJ, K, mapfn,
-                                   degree, tdim, entity, vs, map_type);
+                                   degree, tdim, entity, vs, map_type, ptype);
       transform.insert(transform.end(), t2b.begin(), t2b.end());
     }
 
@@ -468,7 +469,7 @@ doftransforms::compute_entity_transformations(
         4>&,
     std::experimental::mdspan<const float,
                               std::experimental::dextents<std::size_t, 2>>,
-    int, std::size_t, maps::type);
+    int, std::size_t, maps::type, polyset::type);
 
 template std::map<cell::type,
                   std::pair<std::vector<double>, std::array<std::size_t, 3>>>
@@ -484,6 +485,6 @@ doftransforms::compute_entity_transformations(
         4>&,
     std::experimental::mdspan<const double,
                               std::experimental::dextents<std::size_t, 2>>,
-    int, std::size_t, maps::type);
+    int, std::size_t, maps::type, polyset::type);
 /// @endcond
 //-----------------------------------------------------------------------------

cpp/basix/dof-transformations.h

@@ -7,6 +7,7 @@
 #include "cell.h"
 #include "maps.h"
 #include "mdspan.hpp"
+#include "polyset.h"
 #include <array>
 #include <concepts>
 #include <map>
@@ -33,6 +34,7 @@ namespace basix::doftransforms
 /// @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
+/// @param[in] ptype The polyset type used by the element
 /// @return Entity transformations. For each cell, the shape is
 /// (ntransformation, ndofs, ndofs)
 template <std::floating_point T>
@@ -48,6 +50,6 @@ compute_entity_transformations(
     std::experimental::mdspan<const T,
                               std::experimental::dextents<std::size_t, 2>>
         coeffs,
-    int degree, std::size_t vs, maps::type map_type);
+    int degree, std::size_t vs, maps::type map_type, polyset::type ptype);
 
 } // namespace basix::doftransforms

cpp/basix/e-lagrange.cpp

@@ -1434,10 +1434,10 @@ FiniteElement<T> basix::element::create_iso(cell::type celltype, int degree,
                                             lagrange_variant variant,
                                             bool discontinuous)
 {
-  if (celltype != cell::type::interval)
+  if (celltype != cell::type::interval && celltype != cell::type::quadrilateral)
   {
-    throw std::runtime_error(
-        "Can currently only create iso elements on an interval");
+    throw std::runtime_error("Can currently only create iso elements on "
+                             "intervals and quadrilaterals");
   }
 
   if (variant == lagrange_variant::unset)

cpp/basix/finite-element.cpp

@@ -658,7 +658,7 @@ FiniteElement<F>::FiniteElement(
   _entity_transformations = doftransforms::compute_entity_transformations(
       cell_type, x, M,
       mdspan_t<const F, 2>(_coeffs.first.data(), _coeffs.second),
-      highest_degree, value_size, map_type);
+      highest_degree, value_size, map_type, poly_type);
 
   const std::size_t nderivs
       = polyset::nderivs(cell_type, interpolation_nderivs);