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 "precompute.h"
#include <array>
#include <functional>
#include <numeric>
#include <string>
#include <vector>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
#include <xtl/xspan.hpp>
 
namespace basix
{
 
namespace element
{
 
std::tuple<std::array<std::vector<xt::xtensor<double, 2>>, 4>,
           std::array<std::vector<xt::xtensor<double, 3>>, 4>>
make_discontinuous(const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
                   const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
                   int tdim, int value_size);
 
} // namespace element
 
 
class FiniteElement
{
 
public:
  FiniteElement(
      element::family family, cell::type cell_type, int degree,
      const std::vector<std::size_t>& value_shape,
      const xt::xtensor<double, 2>& wcoeffs,
      const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
      const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
      maps::type map_type, bool discontinuous, int highest_complete_degree,
      std::vector<std::tuple<std::vector<FiniteElement>, std::vector<int>>>
          tensor_factors
      = {},
      element::lagrange_variant lvariant = element::lagrange_variant::unset,
      element::dpc_variant dvariant = element::dpc_variant::unset);
 
  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;
 
  xt::xtensor<double, 4> tabulate(int nd, const xt::xarray<double>& x) const;
 
  void tabulate(int nd, const xt::xarray<double>& x,
                xt::xtensor<double, 4>& basis) const;
 
  cell::type cell_type() const;
 
  int degree() const;
 
  const std::vector<int>& 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;
 
  bool discontinuous() const;
 
  bool dof_transformations_are_permutations() const;
 
  bool dof_transformations_are_identity() const;
 
  xt::xtensor<double, 3> push_forward(const xt::xtensor<double, 3>& U,
                                      const xt::xtensor<double, 3>& J,
                                      const xtl::span<const double>& detJ,
                                      const xt::xtensor<double, 3>& K) const;
 
  xt::xtensor<double, 3> pull_back(const xt::xtensor<double, 3>& u,
                                   const xt::xtensor<double, 3>& J,
                                   const xtl::span<const double>& detJ,
                                   const xt::xtensor<double, 3>& 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&) { u.assign(U); };
    case maps::type::L2Piola:
      return [](O& u, const P& U, const Q& J, double detJ, const R& K) {
        maps::l2_piola(u, U, J, detJ, K);
      };
    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<int>>& num_entity_dofs() const;
 
  const std::vector<std::vector<int>>& num_entity_closure_dofs() const;
 
  const std::vector<std::vector<std::vector<int>>>& entity_dofs() const;
 
  const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const;
 
  xt::xtensor<double, 3> base_transformations() const;
 
  std::map<cell::type, xt::xtensor<double, 3>> entity_transformations() const;
 
  void permute_dofs(const xtl::span<std::int32_t>& dofs,
                    std::uint32_t cell_info) const;
 
  void unpermute_dofs(const xtl::span<std::int32_t>& dofs,
                      std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_dof_transformation(const xtl::span<T>& data, int block_size,
                                std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_transpose_dof_transformation(const xtl::span<T>& data,
                                          int block_size,
                                          std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_transpose_dof_transformation(
      const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_dof_transformation(const xtl::span<T>& data,
                                        int block_size,
                                        std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_dof_transformation_to_transpose(const xtl::span<T>& data,
                                             int block_size,
                                             std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_transpose_dof_transformation_to_transpose(
      const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_transpose_dof_transformation_to_transpose(
      const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  template <typename T>
  void apply_inverse_dof_transformation_to_transpose(
      const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
 
  const xt::xtensor<double, 2>& points() const;
 
  const xt::xtensor<double, 2>& interpolation_matrix() const;
 
  const xt::xtensor<double, 2>& dual_matrix() const;
 
  const xt::xtensor<double, 2>& wcoeffs() const;
 
  const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x() const;
 
  const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M() const;
 
  const xt::xtensor<double, 2>& coefficient_matrix() const;
 
  std::array<int, 2> degree_bounds() 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;
 
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
  int _degree;
 
  // Value shape
  std::vector<int> _value_shape;
 
  maps::type _map_type;
 
  // 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$).
  xt::xtensor<double, 2> _coeffs;
 
  // Number of dofs associated with each cell (sub-)entity
  //
  // The dofs of an element are associated with entities of different
  // topological dimension (vertices, edges, faces, cells). The dofs are
  // listed in this order, with vertex dofs first. Each entry is the dof
  // count on the associated entity, as listed by cell::topology.
  std::vector<std::vector<int>> _num_edofs;
 
  // Number of dofs associated with the closure of each cell
  // (sub-)entity
  std::vector<std::vector<int>> _num_e_closure_dofs;
 
  // Dofs associated with each cell (sub-)entity
  std::vector<std::vector<std::vector<int>>> _edofs;
 
  // Dofs associated with each cell (sub-)entity
  std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
 
  // Entity transformations
  std::map<cell::type, xt::xtensor<double, 3>> _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
  xt::xtensor<double, 2> _points;
 
  // Interpolation points on the cell. The shape is (entity_dim, num
  // entities of given dimension, num_points, tdim)
  std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
 
  xt::xtensor<double, 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;
 
  // The entity transformations in precomputed form
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans;
 
  // The transposed entity transformations in precomputed form
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etransT;
 
  // The inverse entity transformations in precomputed form
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans_inv;
 
  // The inverse transpose entity transformations in precomputed form
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans_invT;
 
  // Indicates whether or not this is the discontinuous version of the
  // element
  bool _discontinuous;
 
  // The dual matrix
  xt::xtensor<double, 2> _dual_matrix;
 
  // Polynomial degree bounds
  // [0]: highest degree n such that Lagrange order n is a subspace of
  // this space
  // [1]: highest polynomial degree
  std::array<int, 2> _degree_bounds;
 
  // 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
  xt::xtensor<double, 2> _wcoeffs;
 
  // Interpolation matrices for each entity
  std::array<std::vector<xt::xtensor<double, 3>>, 4> _M;
};
 
FiniteElement create_custom_element(
    cell::type cell_type, int degree,
    const std::vector<std::size_t>& value_shape,
    const xt::xtensor<double, 2>& wcoeffs,
    const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
    const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
    maps::type map_type, bool discontinuous, int highest_complete_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 xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix(_etrans.at(cell::type::interval)[0], data,
                                 dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[1],
                                   data, dofstart, block_size);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[0],
                                   data, dofstart, block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_transpose_dof_transformation(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix(_etransT.at(cell::type::interval)[0], data,
                                 dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[0],
                                   data, dofstart, block_size);
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[1],
                                   data, dofstart, block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_transpose_dof_transformation(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix(_etrans_invT.at(cell::type::interval)[0], data,
                                 dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix(
              _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix(
              _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_dof_transformation(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix(_etrans_inv.at(cell::type::interval)[0], data,
                                 dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix(
              _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix(
              _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_dof_transformation_to_transpose(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix_to_transpose(
            _etrans.at(cell::type::interval)[0], data, dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix_to_transpose(
              _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix_to_transpose(
              _etrans.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_transpose_dof_transformation_to_transpose(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix_to_transpose(
            _etrans_invT.at(cell::type::interval)[0], data, dofstart,
            block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix_to_transpose(
              _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
 
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix_to_transpose(
              _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_transpose_dof_transformation_to_transpose(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix_to_transpose(
            _etransT.at(cell::type::interval)[0], data, dofstart, block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix_to_transpose(
              _etransT.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix_to_transpose(
              _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
 
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::apply_inverse_dof_transformation_to_transpose(
    const xtl::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 * _num_edofs[2].size() : 0;
    int dofstart
        = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
 
    // Transform DOFs on edges
    for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
    {
      // Reverse an edge
      if (cell_info >> (face_start + e) & 1)
        precompute::apply_matrix_to_transpose(
            _etrans_inv.at(cell::type::interval)[0], data, dofstart,
            block_size);
      dofstart += _num_edofs[1][e];
    }
 
    if (_cell_tdim == 3)
    {
      // Permute DOFs on faces
      for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
      {
        // Rotate a face
        for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
          precompute::apply_matrix_to_transpose(
              _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
              block_size);
 
        // Reflect a face
        if (cell_info >> (3 * f) & 1)
          precompute::apply_matrix_to_transpose(
              _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
              block_size);
 
        dofstart += _num_edofs[2][f];
      }
    }
  }
}
//-----------------------------------------------------------------------------
 
} // namespace basix