finite-element.h

// Copyright (c) 2020 Chris Richardson
// FEniCS Project
// SPDX-License-Identifier:    MIT
 
// FIXME: just include everything for now
// Need to define public API
 
#pragma once
 
#include "cell.h"
#include "element-families.h"
#include "lattice.h"
#include "maps.h"
#include "precompute.h"
#include <array>
#include <numeric>
#include <set>
#include <string>
#include <vector>
#include <xtensor/xadapt.hpp>
#include <xtensor/xcomplex.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
#include <xtl/xspan.hpp>
 
namespace basix
{
 
xt::xtensor<double, 3> compute_expansion_coefficients(
    cell::type cell_type, const xt::xtensor<double, 2>& B,
    const std::vector<std::vector<xt::xtensor<double, 3>>>& M,
    const std::vector<std::vector<xt::xtensor<double, 2>>>& x, int degree,
    double kappa_tol = 0.0);
 
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, 3>& coeffs,
                const std::map<cell::type, xt::xtensor<double, 3>>&
                    entity_transformations,
                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 = maps::type::identity);
 
  FiniteElement(const FiniteElement& element) = default;
 
  FiniteElement(FiniteElement&& element) = default;
 
  ~FiniteElement() = default;
 
  FiniteElement& operator=(const FiniteElement& element) = default;
 
  FiniteElement& operator=(FiniteElement&& element) = default;
 
  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_data) const;
 
  cell::type cell_type() const;
 
  int degree() const;
 
  int value_size() const;
 
  const std::vector<int>& value_shape() const;
 
  int dim() const;
 
  element::family family() const;
 
  maps::type mapping_type() const;
 
  bool dof_transformations_are_permutations() const;
 
  bool dof_transformations_are_identity() const;
 
  xt::xtensor<double, 3>
  map_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;
 
  template <typename O, typename P, typename Q, typename S, typename T>
  void map_push_forward_m(const O& U, const P& J, const Q& detJ, const S& K,
                          T&& u) const
  {
    // FIXME: Should U.shape(2) be replaced by the physical value size?
    // Can it differ?
 
    // Loop over points that share J
    for (std::size_t i = 0; i < U.shape(0); ++i)
    {
      auto _J = xt::view(J, i, xt::all(), xt::all());
      auto _K = xt::view(K, i, xt::all(), xt::all());
      auto _U = xt::view(U, i, xt::all(), xt::all());
      auto _u = xt::view(u, i, xt::all(), xt::all());
      maps::apply_map(_u, _U, _J, detJ[i], _K, map_type);
    }
  }
 
  xt::xtensor<double, 3> map_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 S, typename T>
  void map_pull_back_m(const O& u, const P& J, const Q& detJ, const S& K,
                       T&& U) const
  {
    // Loop over points that share K and K
    for (std::size_t i = 0; i < u.shape(0); ++i)
    {
      auto _J = xt::view(J, i, xt::all(), xt::all());
      auto _K = xt::view(K, i, xt::all(), xt::all());
      auto _u = xt::view(u, i, xt::all(), xt::all());
      auto _U = xt::view(U, i, xt::all(), xt::all());
      maps::apply_map(_U, _u, _K, 1.0 / detJ[i], _J, map_type);
    }
  }
 
  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::set<int>>>& entity_dofs() const;
 
  const std::vector<std::vector<std::set<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;
 
  int num_points() const;
 
  const xt::xtensor<double, 2>& interpolation_matrix() const;
 
  template <typename T>
  void interpolate(const xtl::span<T>& coefficients,
                   const xtl::span<const T>& data, const int block_size) const;
 
  maps::type map_type;
 
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;
 
  // 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$. i.e.,
  // _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::set<int>>> _edofs;
 
  // Dofs associated with each cell (sub-)entity
  std::vector<std::vector<std::set<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;
 
  std::array<std::vector<xt::xtensor<double, 3>>, 4> _matM_new;
 
  bool _dof_transformations_are_permutations;
 
  bool _dof_transformations_are_identity;
 
  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
 
  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
 
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans;
 
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etransT;
 
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans_inv;
 
  std::map<cell::type,
           std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
                                  xt::xtensor<double, 2>>>>
      _etrans_invT;
};
 
FiniteElement
create_element(element::family family, cell::type cell, int degree,
               lattice::type lattice_type);
 
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];
      }
    }
  }
}
//-----------------------------------------------------------------------------
template <typename T>
void FiniteElement::interpolate(const xtl::span<T>& coefficients,
                                const xtl::span<const T>& data,
                                const int block_size) const
{
  if (block_size != 1)
  {
    throw std::runtime_error(
        "Interpolation of blocked data not implemented yet.");
  }
 
  const std::size_t rows = dim();
 
  // Compute coefficients = Pi * x (matrix-vector multiply)
  const xt::xtensor<double, 2>& Pi = interpolation_matrix();
  assert(Pi.size() % rows == 0);
  const std::size_t cols = Pi.size() / rows;
  for (std::size_t i = 0; i < rows; ++i)
  {
    // Can be replaced with std::transform_reduce once GCC 8 series dies.
    // Dot product between row i of the matrix and 'data'
    coefficients[i] = std::inner_product(std::next(Pi.data(), i * cols),
                                         std::next(Pi.data(), i * cols + cols),
                                         data.data(), T(0.0));
  }
}
//-----------------------------------------------------------------------------
 
} // namespace basix