finite-element.h

// Copyright (c) 2020 Chris Richardson
// FEniCS Project
// SPDX-License-Identifier:    MIT
 
#pragma once
 
#include "cell.h"
#include "element-families.h"
#include "maps.h"
#include "mdspan.hpp"
#include "precompute.h"
#include "sobolev-spaces.h"
#include <array>
#include <functional>
#include <map>
#include <numeric>
#include <span>
#include <string>
#include <tuple>
#include <vector>
 
namespace basix
{
 
namespace impl
{
using mdspan2_t
    = std::experimental::mdspan<double,
                                std::experimental::dextents<std::size_t, 2>>;
using mdspan4_t
    = std::experimental::mdspan<double,
                                std::experimental::dextents<std::size_t, 4>>;
using cmdspan2_t
    = std::experimental::mdspan<const double,
                                std::experimental::dextents<std::size_t, 2>>;
using cmdspan3_t
    = std::experimental::mdspan<const double,
                                std::experimental::dextents<std::size_t, 3>>;
using cmdspan4_t
    = std::experimental::mdspan<const double,
                                std::experimental::dextents<std::size_t, 4>>;
 
using mdarray2_t
    = std::experimental::mdarray<double,
                                 std::experimental::dextents<std::size_t, 2>>;
using mdarray4_t
    = std::experimental::mdarray<double,
                                 std::experimental::dextents<std::size_t, 4>>;
 
inline std::array<std::vector<cmdspan2_t>, 4>
to_mdspan(std::array<std::vector<mdarray2_t>, 4>& x)
{
  std::array<std::vector<cmdspan2_t>, 4> x1;
  for (std::size_t i = 0; i < x.size(); ++i)
    for (std::size_t j = 0; j < x[i].size(); ++j)
      x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
 
  return x1;
}
 
inline std::array<std::vector<cmdspan4_t>, 4>
to_mdspan(std::array<std::vector<mdarray4_t>, 4>& M)
{
  std::array<std::vector<cmdspan4_t>, 4> M1;
  for (std::size_t i = 0; i < M.size(); ++i)
    for (std::size_t j = 0; j < M[i].size(); ++j)
      M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
 
  return M1;
}
 
inline std::array<std::vector<cmdspan2_t>, 4>
to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& x,
          const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
{
  std::array<std::vector<cmdspan2_t>, 4> x1;
  for (std::size_t i = 0; i < x.size(); ++i)
    for (std::size_t j = 0; j < x[i].size(); ++j)
      x1[i].push_back(cmdspan2_t(x[i][j].data(), shape[i][j]));
 
  return x1;
}
 
inline std::array<std::vector<cmdspan4_t>, 4>
to_mdspan(const std::array<std::vector<std::vector<double>>, 4>& M,
          const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
{
  std::array<std::vector<cmdspan4_t>, 4> M1;
  for (std::size_t i = 0; i < M.size(); ++i)
    for (std::size_t j = 0; j < M[i].size(); ++j)
      M1[i].push_back(cmdspan4_t(M[i][j].data(), shape[i][j]));
 
  return M1;
}
 
} // namespace impl
 
namespace element
{
using cmdspan2_t = impl::cmdspan2_t;
using cmdspan4_t = impl::cmdspan4_t;
 
std::tuple<std::array<std::vector<std::vector<double>>, 4>,
           std::array<std::vector<std::array<std::size_t, 2>>, 4>,
           std::array<std::vector<std::vector<double>>, 4>,
           std::array<std::vector<std::array<std::size_t, 4>>, 4>>
make_discontinuous(const std::array<std::vector<cmdspan2_t>, 4>& x,
                   const std::array<std::vector<cmdspan4_t>, 4>& M,
                   std::size_t tdim, std::size_t value_size);
 
} // namespace element
 
class FiniteElement
{
  using cmdspan2_t
      = std::experimental::mdspan<const double,
                                  std::experimental::dextents<std::size_t, 2>>;
  using cmdspan4_t
      = std::experimental::mdspan<const double,
                                  std::experimental::dextents<std::size_t, 4>>;
 
public:
  FiniteElement(
      element::family family, cell::type cell_type, int degree,
      const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
      const std::array<std::vector<cmdspan2_t>, 4>& x,
      const std::array<std::vector<cmdspan4_t>, 4>& M,
      int interpolation_nderivs, maps::type map_type,
      sobolev::space sobolev_space, bool discontinuous,
      int highest_complete_degree, int highest_degree,
      element::lagrange_variant lvariant, element::dpc_variant dvariant,
      std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
          tensor_factors
      = {});
 
  FiniteElement(
      element::family family, cell::type cell_type, int degree,
      const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
      const std::array<std::vector<cmdspan2_t>, 4>& x,
      const std::array<std::vector<cmdspan4_t>, 4>& M,
      int interpolation_nderivs, maps::type map_type,
      sobolev::space sobolev_space, bool discontinuous,
      int highest_complete_degree, int highest_degree,
      element::lagrange_variant lvariant,
      std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
          tensor_factors
      = {});
 
  FiniteElement(
      element::family family, cell::type cell_type, int degree,
      const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
      const std::array<std::vector<cmdspan2_t>, 4>& x,
      const std::array<std::vector<cmdspan4_t>, 4>& M,
      int interpolation_nderivs, maps::type map_type,
      sobolev::space sobolev_space, bool discontinuous,
      int highest_complete_degree, int highest_degree,
      element::dpc_variant dvariant,
      std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
          tensor_factors
      = {});
 
  FiniteElement(
      element::family family, cell::type cell_type, int degree,
      const std::vector<std::size_t>& value_shape, const cmdspan2_t& wcoeffs,
      const std::array<std::vector<cmdspan2_t>, 4>& x,
      const std::array<std::vector<cmdspan4_t>, 4>& M,
      int interpolation_nderivs, maps::type map_type,
      sobolev::space sobolev_space, bool discontinuous,
      int highest_complete_degree, int highest_degree,
      std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
          tensor_factors
      = {});
 
  FiniteElement(const FiniteElement& element) = default;
 
  FiniteElement(FiniteElement&& element) = default;
 
  ~FiniteElement() = default;
 
  FiniteElement& operator=(const FiniteElement& element) = default;
 
  FiniteElement& operator=(FiniteElement&& element) = default;
 
  bool operator==(const FiniteElement& e) const;
 
  std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
                                            std::size_t num_points) const;
 
  std::pair<std::vector<double>, std::array<std::size_t, 4>>
  tabulate(int nd, impl::cmdspan2_t x) const;
 
  std::pair<std::vector<double>, std::array<std::size_t, 4>>
  tabulate(int nd, const std::span<const double>& x,
           std::array<std::size_t, 2> shape) const;
 
  void tabulate(int nd, impl::cmdspan2_t x, impl::mdspan4_t basis) const;
 
  void tabulate(int nd, const std::span<const double>& x,
                std::array<std::size_t, 2> xshape,
                const std::span<double>& basis) const;
 
  cell::type cell_type() const;
 
  int degree() const;
 
  int highest_degree() const;
 
  int highest_complete_degree() const;
 
  const std::vector<std::size_t>& value_shape() const;
 
  int dim() const;
 
  element::family family() const;
 
  element::lagrange_variant lagrange_variant() const;
 
  element::dpc_variant dpc_variant() const;
 
  maps::type map_type() const;
 
  sobolev::space sobolev_space() const;
 
  bool discontinuous() const;
 
  bool dof_transformations_are_permutations() const;
 
  bool dof_transformations_are_identity() const;
 
  std::pair<std::vector<double>, std::array<std::size_t, 3>>
  push_forward(impl::cmdspan3_t U, impl::cmdspan3_t J,
               std::span<const double> detJ, impl::cmdspan3_t K) const;
 
  std::pair<std::vector<double>, std::array<std::size_t, 3>>
  pull_back(impl::cmdspan3_t u, impl::cmdspan3_t J,
            std::span<const double> detJ, impl::cmdspan3_t K) const;
 
  template <typename O, typename P, typename Q, typename R>
  std::function<void(O&, const P&, const Q&, double, const R&)> map_fn() const
  {
    switch (_map_type)
    {
    case maps::type::identity:
      return [](O& u, const P& U, const Q&, double, const R&)
      {
        assert(U.extent(0) == u.extent(0));
        assert(U.extent(1) == u.extent(1));
        for (std::size_t i = 0; i < U.extent(0); ++i)
          for (std::size_t j = 0; j < U.extent(1); ++j)
            u(i, j) = U(i, j);
      };
    case maps::type::covariantPiola:
      return [](O& u, const P& U, const Q& J, double detJ, const R& K)
      { maps::covariant_piola(u, U, J, detJ, K); };
    case maps::type::contravariantPiola:
      return [](O& u, const P& U, const Q& J, double detJ, const R& K)
      { maps::contravariant_piola(u, U, J, detJ, K); };
    case maps::type::doubleCovariantPiola:
      return [](O& u, const P& U, const Q& J, double detJ, const R& K)
      { maps::double_covariant_piola(u, U, J, detJ, K); };
    case maps::type::doubleContravariantPiola:
      return [](O& u, const P& U, const Q& J, double detJ, const R& K)
      { maps::double_contravariant_piola(u, U, J, detJ, K); };
    default:
      throw std::runtime_error("Map not implemented");
    }
  }
 
  const std::vector<std::vector<std::vector<int>>>& entity_dofs() const;
 
  const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const;
 
  std::pair<std::vector<double>, std::array<std::size_t, 3>>
  base_transformations() const;
 
  std::map<cell::type,
           std::pair<std::vector<double>, std::array<std::size_t, 3>>>
  entity_transformations() const;
 
  void permute_dofs(const std::span<std::int32_t>& dofs,
                    std::uint32_t cell_info) const;
 
  void unpermute_dofs(const std::span<std::int32_t>& dofs,
                      std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_dof_transformation(const std::span<T>& data, int block_size,
                                std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_transpose_dof_transformation(const std::span<T>& data,
                                          int block_size,
                                          std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_transpose_dof_transformation(
      const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_dof_transformation(const std::span<T>& data,
                                        int block_size,
                                        std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_dof_transformation_to_transpose(const std::span<T>& data,
                                             int block_size,
                                             std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_transpose_dof_transformation_to_transpose(
      const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_transpose_dof_transformation_to_transpose(
      const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_dof_transformation_to_transpose(
      const std::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
  points() const;
 
  const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
  interpolation_matrix() const;
 
  const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
  dual_matrix() const;
 
  const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
  wcoeffs() const;
 
  const std::array<
      std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
      4>&
  x() const;
 
  const std::array<
      std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>,
      4>&
  M() const;
 
  const std::pair<std::vector<double>, std::array<std::size_t, 2>>&
  coefficient_matrix() const;
 
  bool has_tensor_product_factorisation() const;
 
  std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
  get_tensor_product_representation() const;
 
  bool interpolation_is_identity() const;
 
  int interpolation_nderivs() const;
 
private:
  // Cell type
  cell::type _cell_type;
 
  // Topological dimension of the cell
  std::size_t _cell_tdim;
 
  // Topological dimension of the cell
  std::vector<std::vector<cell::type>> _cell_subentity_types;
 
  // Finite element family
  element::family _family;
 
  // Lagrange variant
  element::lagrange_variant _lagrange_variant;
 
  // Lagrange variant
  element::dpc_variant _dpc_variant;
 
  // Degree that was input when creating the element
  int _degree;
 
  // Degree
  int _interpolation_nderivs;
 
  // Highest degree polynomial in element's polyset
  int _highest_degree;
 
  // Highest degree space that is a subspace of element's polyset
  int _highest_complete_degree;
 
  // Value shape
  std::vector<std::size_t> _value_shape;
 
  maps::type _map_type;
 
  sobolev::space _sobolev_space;
 
  // Shape function coefficient of expansion sets on cell. If shape
  // function is given by @f$\psi_i = \sum_{k} \phi_{k}
  // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
  // _coeffs.row(i) are the expansion coefficients for shape function i
  // (@f$\psi_{i}@f$).
  std::pair<std::vector<double>, std::array<std::size_t, 2>> _coeffs;
 
  // Dofs associated with each cell (sub-)entity
  std::vector<std::vector<std::vector<int>>> _edofs;
 
  // Dofs associated with the closdure of each cell (sub-)entity
  std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
 
  using array2_t = std::pair<std::vector<double>, std::array<std::size_t, 2>>;
  using array3_t = std::pair<std::vector<double>, std::array<std::size_t, 3>>;
 
  // Entity transformations
  std::map<cell::type, array3_t> _entity_transformations;
 
  // Set of points used for point evaluation
  // Experimental - currently used for an implementation of
  // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
  // away. For non-Lagrange elements, these points will be used in combination
  // with _interpolation_matrix to perform interpolation
  std::pair<std::vector<double>, std::array<std::size_t, 2>> _points;
 
  // Interpolation points on the cell. The shape is (entity_dim, num
  // entities of given dimension, num_points, tdim)
  std::array<
      std::vector<std::pair<std::vector<double>, std::array<std::size_t, 2>>>,
      4>
      _x;
 
  std::pair<std::vector<double>, std::array<std::size_t, 2>> _matM;
 
  // Indicates whether or not the DOF transformations are all
  // permutations
  bool _dof_transformations_are_permutations;
 
  // Indicates whether or not the DOF transformations are all identity
  bool _dof_transformations_are_identity;
 
  // The entity permutations (factorised). This will only be set if
  // _dof_transformations_are_permutations is True and
  // _dof_transformations_are_identity is False
  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
 
  // The reverse entity permutations (factorised). This will only be set
  // if _dof_transformations_are_permutations is True and
  // _dof_transformations_are_identity is False
  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
 
  using trans_data_t
      = std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
 
  // The entity transformations in precomputed form
  std::map<cell::type, trans_data_t> _etrans;
 
  // The transposed entity transformations in precomputed form
  std::map<cell::type, trans_data_t> _etransT;
 
  // The inverse entity transformations in precomputed form
  std::map<cell::type, trans_data_t> _etrans_inv;
 
  // The inverse transpose entity transformations in precomputed form
  std::map<cell::type, trans_data_t> _etrans_invT;
 
  // Indicates whether or not this is the discontinuous version of the
  // element
  bool _discontinuous;
 
  // The dual matrix
  std::pair<std::vector<double>, std::array<std::size_t, 2>> _dual_matrix;
 
  // Tensor product representation
  // Entries of tuple are (list of elements on an interval, permutation
  // of DOF numbers)
  // @todo: For vector-valued elements, a tensor product type and a
  // scaling factor may additionally be needed.
  std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
      _tensor_factors;
 
  // Is the interpolation matrix an identity?
  bool _interpolation_is_identity;
 
  // The coefficients that define the polynomial set in terms of the
  // orthonormal polynomials
  std::pair<std::vector<double>, std::array<std::size_t, 2>> _wcoeffs;
 
  // Interpolation matrices for each entity
  using array4_t
      = std::vector<std::pair<std::vector<double>, std::array<std::size_t, 4>>>;
  std::array<array4_t, 4> _M;
  // std::array<
  //     std::vector<std::pair<std::vector<double>, std::array<std::size_t,
  //     4>>>, 4> _M;
};
 
FiniteElement
create_custom_element(cell::type cell_type,
                      const std::vector<std::size_t>& value_shape,
                      const impl::cmdspan2_t& wcoeffs,
                      const std::array<std::vector<impl::cmdspan2_t>, 4>& x,
                      const std::array<std::vector<impl::cmdspan4_t>, 4>& M,
                      int interpolation_nderivs, maps::type map_type,
                      sobolev::space sobolev_space, bool discontinuous,
                      int highest_complete_degree, int highest_degree);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::lagrange_variant lvariant,
                             bool discontinuous);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::lagrange_variant lvariant,
                             element::dpc_variant dvariant, bool discontinuous);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::dpc_variant dvariant,
                             bool discontinuous);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, bool discontinuous);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::lagrange_variant lvariant);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::lagrange_variant lvariant,
                             element::dpc_variant dvariant);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree, element::dpc_variant dvariant);
 
FiniteElement create_element(element::family family, cell::type cell,
                             int degree);
 
std::string version();
 
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_dof_transformation(const std::span<T>& data,
                                             int block_size,
                                             std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans.at(_cell_subentity_types[2][f]);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          const auto& m = trans[1];
          const auto& v_size_t = std::get<0>(m);
          const auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          const auto& m = trans[0];
          const auto& v_size_t = std::get<0>(m);
          const auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_transpose_dof_transformation(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etransT.at(_cell_subentity_types[2][f]);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_transpose_dof_transformation(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
      {
        auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
        precompute::apply_matrix(std::span(v_size_t),
                                 cmdspan2_t(matrix.first.data(), matrix.second),
                                 data, dofstart, block_size);
      }
      dofstart += _edofs[1][e].size();
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          const auto& m = trans[0];
          const auto& v_size_t = std::get<0>(m);
          const auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_dof_transformation(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_dof_transformation_to_transpose(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etrans.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans.at(_cell_subentity_types[2][f]);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_transpose_dof_transformation_to_transpose(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etrans_invT.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans_invT.at(_cell_subentity_types[2][f]);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_transpose_dof_transformation_to_transpose(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etransT.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etransT.at(_cell_subentity_types[2][f]);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_dof_transformation_to_transpose(
    const std::span<T>& data, int block_size, std::uint32_t cell_info) const
{
  if (_dof_transformations_are_identity)
    return;
 
  if (_cell_tdim >= 2)
  {
    // This assumes 3 bits are used per face. This will need updating if
    // 3D cells with faces with more than 4 sides are implemented
    int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
    int dofstart = 0;
    for (auto& edofs0 : _edofs[0])
      dofstart += edofs0.size();
 
    // Transform DOFs on edges
    {
      auto& [v_size_t, matrix] = _etrans_inv.at(cell::type::interval)[0];
      for (std::size_t e = 0; e < _edofs[1].size(); ++e)
      {
        // Reverse an edge
        if (cell_info >> (face_start + e) & 1)
        {
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
        dofstart += _edofs[1][e].size();
      }
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _edofs[2].size(); ++f)
      {
        auto& trans = _etrans_inv.at(_cell_subentity_types[2][f]);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
        {
          auto& m = trans[0];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
        {
          auto& m = trans[1];
          auto& v_size_t = std::get<0>(m);
          auto& matrix = std::get<1>(m);
          precompute::apply_matrix_to_transpose(
              std::span(v_size_t),
              cmdspan2_t(matrix.first.data(), matrix.second), data, dofstart,
              block_size);
        }
 
        dofstart += _edofs[2][f].size();
      }
    }
  }
}
//-----------------------------------------------------------------------------
 
} // namespace basix