FiniteElement.h

// Copyright (C) 2008-2020 Anders Logg and Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <basix/finite-element.h>
#include <dolfinx/common/types.h>
#include <dolfinx/mesh/cell_types.h>
#include <functional>
#include <memory>
#include <numeric>
#include <vector>
#include <xtl/xspan.hpp>
 
struct ufc_finite_element;
 
namespace dolfinx::fem
{
class FiniteElement
{
public:
  explicit FiniteElement(const ufc_finite_element& ufc_element);
 
  FiniteElement(const FiniteElement& element) = delete;
 
  FiniteElement(FiniteElement&& element) = default;
 
  virtual ~FiniteElement() = default;
 
  FiniteElement& operator=(const FiniteElement& element) = delete;
 
  FiniteElement& operator=(FiniteElement&& element) = default;
 
  std::string signature() const noexcept;
 
  mesh::CellType cell_shape() const noexcept;
 
  int space_dimension() const noexcept;
 
  int block_size() const noexcept;
 
  int value_size() const noexcept;
 
  int reference_value_size() const noexcept;
 
  int value_rank() const noexcept;
 
  int value_dimension(int i) const;
 
  std::string family() const noexcept;
 
  // reference_values[num_points][num_dofs][reference_value_size]
  void evaluate_reference_basis(xt::xtensor<double, 3>& values,
                                const xt::xtensor<double, 2>& X) const;
 
  // reference_value_derivatives[num_points][num_dofs][reference_value_size][num_derivatives]
  // void
  // evaluate_reference_basis_derivatives(std::vector<double>& reference_values,
  //                                      int order,
  //                                      const xt::xtensor<double, 2>& X)
  //                                      const;
 
  void transform_reference_basis(xt::xtensor<double, 3>& values,
                                 const xt::xtensor<double, 3>& reference_values,
                                 const xt::xtensor<double, 3>& J,
                                 const xtl::span<const double>& detJ,
                                 const xt::xtensor<double, 3>& K) const;
 
  void transform_reference_basis_derivatives(
      std::vector<double>& values, std::size_t order,
      const std::vector<double>& reference_values, const array2d<double>& X,
      const std::vector<double>& J, const xtl::span<const double>& detJ,
      const std::vector<double>& K) const;
 
  int num_sub_elements() const noexcept;
 
  const std::vector<std::shared_ptr<const FiniteElement>>&
  sub_elements() const noexcept;
 
  std::size_t hash() const noexcept;
 
  std::shared_ptr<const FiniteElement>
  extract_sub_element(const std::vector<int>& component) const;
 
  bool interpolation_ident() const noexcept;
 
  const xt::xtensor<double, 2>& interpolation_points() const;
 
  template <typename T>
  constexpr void interpolate(const xt::xtensor<T, 2>& values,
                             xtl::span<T> dofs) const
  {
    if (!_element)
    {
      throw std::runtime_error("No underlying element for interpolation. "
                               "Cannot interpolate mixed elements directly.");
    }
 
    const std::size_t rows = _space_dim / _bs;
    assert(_space_dim % _bs == 0);
    assert(dofs.size() == rows);
 
    // Compute dofs = Pi * x (matrix-vector multiply)
    const xt::xtensor<double, 2>& Pi = _element->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 'values'
      dofs[i] = std::inner_product(std::next(Pi.data(), i * cols),
                                   std::next(Pi.data(), i * cols + cols),
                                   values.data(), T(0.0));
    }
  }
 
  bool needs_permutation_data() const noexcept;
 
  template <typename T>
  void apply_dof_transformation(xtl::span<T> data,
                                std::uint32_t cell_permutation,
                                int block_size) const
  {
    assert(_element);
    _element->apply_dof_transformation(data, block_size, cell_permutation);
  }
 
  template <typename T>
  void apply_inverse_transpose_dof_transformation(
      xtl::span<T> data, std::uint32_t cell_permutation, int block_size) const
  {
    assert(_element);
    _element->apply_inverse_transpose_dof_transformation(data, block_size,
                                                         cell_permutation);
  }
 
  template <typename T>
  void
  map_pull_back(const xt::xtensor<T, 3>& u, const xt::xtensor<double, 3>& J,
                const xtl::span<const double>& detJ,
                const xt::xtensor<double, 3>& K, xt::xtensor<T, 3>& U) const
  {
    assert(_element);
    _element->map_pull_back_m(u, J, detJ, K, U);
  }
 
private:
  std::string _signature, _family;
 
  mesh::CellType _cell_shape;
 
  int _tdim, _space_dim, _value_size, _reference_value_size;
 
  // List of sub-elements (if any)
  std::vector<std::shared_ptr<const FiniteElement>> _sub_elements;
 
  // Simple hash of the signature string
  std::size_t _hash;
 
  // Dimension of each value space
  std::vector<int> _value_dimension;
 
  // Block size for VectorElements and TensorElements. This gives the
  // number of DOFs colocated at each point.
  int _bs;
 
  // True if element needs dof permutation
  bool _needs_permutation_data;
 
  // Basix Element (nullptr for mixed elements)
  std::unique_ptr<basix::FiniteElement> _element;
};
} // namespace dolfinx::fem