CoordinateElement.h

// Copyright (C) 2018-2020 Garth N. Wells and Chris N. Richardson
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include "ElementDofLayout.h"
#include <basix/element-families.h>
#include <cstdint>
#include <dolfinx/common/math.h>
#include <dolfinx/mesh/cell_types.h>
#include <memory>
#include <xtensor/xtensor.hpp>
#include <xtl/xspan.hpp>
 
namespace basix
{
class FiniteElement;
}
 
namespace dolfinx::fem
{
 
class CoordinateElement
{
public:
  explicit CoordinateElement(
      std::shared_ptr<const basix::FiniteElement> element);
 
  CoordinateElement(mesh::CellType celltype, int degree,
                    basix::element::lagrange_variant type
                    = basix::element::lagrange_variant::equispaced);
 
  virtual ~CoordinateElement() = default;
 
  mesh::CellType cell_shape() const;
 
  int degree() const;
 
  int dim() const;
 
  basix::element::lagrange_variant variant() 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::xtensor<double, 2>& X) const;
 
  void tabulate(int nd, const xt::xtensor<double, 2>& X,
                xt::xtensor<double, 4>& basis) const;
 
  template <typename U, typename V, typename W>
  static void compute_jacobian(const U& dphi, const V& cell_geometry, W&& J)
  {
    math::dot(cell_geometry, dphi, J, true);
  }
 
  template <typename U, typename V>
  static void compute_jacobian_inverse(const U& J, V&& K)
  {
    const int gdim = J.shape(1);
    const int tdim = K.shape(1);
    if (gdim == tdim)
      math::inv(J, K);
    else
      math::pinv(J, K);
  }
 
  template <typename U>
  static double compute_jacobian_determinant(const U& J)
  {
    if (J.shape(0) == J.shape(1))
      return math::det(J);
    else
    {
      using T = typename U::value_type;
      auto B = xt::transpose(J);
      xt::xtensor<T, 2> BA = xt::zeros<T>({B.shape(0), J.shape(1)});
      math::dot(B, J, BA);
      return std::sqrt(math::det(BA));
    }
  }
 
  ElementDofLayout create_dof_layout() const;
 
  static void push_forward(xt::xtensor<double, 2>& x,
                           const xt::xtensor<double, 2>& cell_geometry,
                           const xt::xtensor<double, 2>& phi);
 
  static std::array<double, 3> x0(const xt::xtensor<double, 2>& cell_geometry);
 
  static void pull_back_affine(xt::xtensor<double, 2>& X,
                               const xt::xtensor<double, 2>& K,
                               const std::array<double, 3>& x0,
                               const xt::xtensor<double, 2>& x);
 
  void pull_back_nonaffine(xt::xtensor<double, 2>& X,
                           const xt::xtensor<double, 2>& x,
                           const xt::xtensor<double, 2>& cell_geometry,
                           double tol = 1.0e-8, int maxit = 10) const;
 
  void permute_dofs(const xtl::span<std::int32_t>& dofs,
                    std::uint32_t cell_perm) const;
 
  void unpermute_dofs(const xtl::span<std::int32_t>& dofs,
                      std::uint32_t cell_perm) const;
 
  bool needs_dof_permutations() const;
 
  bool is_affine() const noexcept { return _is_affine; }
 
private:
  // Flag denoting affine map
  bool _is_affine;
 
  // Basix Element
  std::shared_ptr<const basix::FiniteElement> _element;
};
} // namespace dolfinx::fem