Basix
finite-element.h
1 // Copyright (c) 2020 Chris Richardson
2 // FEniCS Project
3 // SPDX-License-Identifier: MIT
4 
5 // FIXME: just include everything for now
6 // Need to define public API
7 
8 #pragma once
9 
10 #include "cell.h"
11 #include "element-families.h"
12 #include "maps.h"
13 #include "span.hpp"
14 #include <array>
15 #include <string>
16 #include <vector>
17 #include <xtensor/xadapt.hpp>
18 #include <xtensor/xcomplex.hpp>
19 #include <xtensor/xtensor.hpp>
20 #include <xtensor/xview.hpp>
21 
23 namespace basix
24 {
25 
181 xt::xtensor<double, 3> compute_expansion_coefficients(
182  cell::type cell_type, const xt::xtensor<double, 2>& B,
183  const std::vector<std::vector<xt::xtensor<double, 3>>>& M,
184  const std::vector<std::vector<xt::xtensor<double, 2>>>& x, int degree,
185  double kappa_tol = 0.0);
186 
191 {
192 
193 public:
209  element::family family, cell::type cell_type, int degree,
210  const std::vector<std::size_t>& value_shape,
211  const xt::xtensor<double, 3>& coeffs,
212  const std::vector<xt::xtensor<double, 2>>& entity_transformations,
213  const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
214  const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
215  maps::type map_type = maps::type::identity);
216 
218  FiniteElement(const FiniteElement& element) = default;
219 
221  FiniteElement(FiniteElement&& element) = default;
222 
224  ~FiniteElement() = default;
225 
227  FiniteElement& operator=(const FiniteElement& element) = default;
228 
230  FiniteElement& operator=(FiniteElement&& element) = default;
231 
249  xt::xtensor<double, 4> tabulate(int nd, const xt::xarray<double>& x) const;
250 
255  void tabulate(int nd, const xt::xarray<double>& x, double* basis_data) const;
256 
259  cell::type cell_type() const;
260 
263  int degree() const;
264 
268  int value_size() const;
269 
272  const std::vector<int>& value_shape() const;
273 
277  int dim() const;
278 
281  element::family family() const;
282 
285  maps::type mapping_type() const;
286 
290 
294 
301  xt::xtensor<double, 3>
302  map_push_forward(const xt::xtensor<double, 3>& U,
303  const xt::xtensor<double, 3>& J,
304  const tcb::span<const double>& detJ,
305  const xt::xtensor<double, 3>& K) const;
306 
313  template <typename T, typename E>
314  void map_push_forward_m(const xt::xtensor<T, 3>& U,
315  const xt::xtensor<double, 3>& J,
316  const tcb::span<const double>& detJ,
317  const xt::xtensor<double, 3>& K, E&& u) const;
318 
325  xt::xtensor<double, 3> map_pull_back(const xt::xtensor<double, 3>& u,
326  const xt::xtensor<double, 3>& J,
327  const tcb::span<const double>& detJ,
328  const xt::xtensor<double, 3>& K) const;
329 
337  template <typename T, typename E>
338  void map_pull_back_m(const xt::xtensor<T, 3>& u,
339  const xt::xtensor<double, 3>& J,
340  const tcb::span<const double>& detJ,
341  const xt::xtensor<double, 3>& K, E&& U) const;
342 
350  const std::vector<std::vector<int>>& entity_dofs() const;
351 
430  xt::xtensor<double, 3> base_transformations() const;
431 
435  void permute_dofs(tcb::span<std::int32_t>& dofs,
436  std::uint32_t cell_info) const;
437 
441  void unpermute_dofs(tcb::span<std::int32_t>& dofs,
442  std::uint32_t cell_info) const;
443 
448  const xt::xtensor<double, 2>& points() const;
449 
451  int num_points() const;
452 
459  const xt::xtensor<double, 2>& interpolation_matrix() const;
460 
463 
464 private:
465  // Cell type
466  cell::type _cell_type;
467 
468  // Topological dimension of the cell
469  std::size_t _cell_tdim;
470 
471  // Finite element family
472  element::family _family;
473 
474  // Degree
475  int _degree;
476 
477  // Value shape
478  std::vector<int> _value_shape;
479 
481  maps::type _map_type;
482 
483  // Shape function coefficient of expansion sets on cell. If shape
484  // function is given by @f$\psi_i = \sum_{k} \phi_{k}
485  // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. i.e.,
486  // _coeffs.row(i) are the expansion coefficients for shape function i
487  // (@f$\psi_{i}@f$).
488  xt::xtensor<double, 2> _coeffs;
489 
490  // Number of cell subentities of each dimension
491  std::vector<int> _cell_sub_entity_count;
492 
493  // Number of dofs associated each subentity
494  //
495  // The dofs of an element are associated with entities of different
496  // topological dimension (vertices, edges, faces, cells). The dofs are
497  // listed in this order, with vertex dofs first. Each entry is the dof
498  // count on the associated entity, as listed by cell::topology.
499  std::vector<std::vector<int>> _entity_dofs;
500 
501  // Entity transformations
502  std::vector<xt::xtensor<double, 2>> _entity_transformations;
503 
504  // Set of points used for point evaluation
505  // Experimental - currently used for an implementation of
506  // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
507  // away. For non-Lagrange elements, these points will be used in combination
508  // with _interpolation_matrix to perform interpolation
509  xt::xtensor<double, 2> _points;
510 
511  // Interpolation points on the cell. The shape is (entity_dim, num
512  // entities of given dimension, num_points, tdim)
513  std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
514 
516  xt::xtensor<double, 2> _matM;
517 
519  std::array<std::vector<xt::xtensor<double, 3>>, 4> _matM_new;
520 
522  bool _dof_transformations_are_permutations;
523 
525  bool _dof_transformations_are_identity;
526 
529  std::vector<std::vector<int>> _entity_permutations;
530 
533  std::vector<std::vector<int>> _reverse_entity_permutations;
534 };
535 
537 FiniteElement create_element(std::string family, std::string cell, int degree);
538 
540 FiniteElement create_element(element::family family, cell::type cell,
541  int degree);
542 
545 std::string version();
546 
547 //-----------------------------------------------------------------------------
548 template <typename T, typename E>
549 void FiniteElement::map_push_forward_m(const xt::xtensor<T, 3>& U,
550  const xt::xtensor<double, 3>& J,
551  const tcb::span<const double>& detJ,
552  const xt::xtensor<double, 3>& K,
553  E&& u) const
554 {
555  // FIXME: Should U.shape(2) be replaced by the physical value size?
556  // Can it differ?
557  // std::array<std::size_t, 3> s = {U.shape(0), U.shape(1), U.shape(2)};
558  // auto _u = xt::adapt(u, s[0] * s[1] * s[2], xt::no_ownership(), s);
559 
560  // Loop over each point
561  for (std::size_t p = 0; p < U.shape(0); ++p)
562  {
563  auto J_p = xt::view(J, p, xt::all(), xt::all());
564  auto K_p = xt::view(K, p, xt::all(), xt::all());
565 
566  // Loop over values at each point
567  for (std::size_t i = 0; i < U.shape(1); ++i)
568  {
569  auto U_data = xt::view(U, p, i, xt::all());
570  auto u_data = xt::view(u, p, i, xt::all());
571  maps::apply_map(u_data, U_data, J_p, detJ[p], K_p, map_type);
572  }
573  }
574 }
575 //-----------------------------------------------------------------------------
576 template <typename T, typename E>
577 void FiniteElement::map_pull_back_m(const xt::xtensor<T, 3>& u,
578  const xt::xtensor<double, 3>& J,
579  const tcb::span<const double>& detJ,
580  const xt::xtensor<double, 3>& K,
581  E&& U) const
582 {
583  // Loop over each point
584  for (std::size_t p = 0; p < u.shape(0); ++p)
585  {
586  auto J_p = xt::view(J, p, xt::all(), xt::all());
587  auto K_p = xt::view(K, p, xt::all(), xt::all());
588 
589  // Loop over each item at point to be transformed
590  for (std::size_t i = 0; i < u.shape(1); ++i)
591  {
592  auto u_data = xt::view(u, p, i, xt::all());
593  auto U_data = xt::view(U, p, i, xt::all());
594  maps::apply_map(U_data, u_data, K_p, 1.0 / detJ[p], J_p, map_type);
595  }
596  }
597 }
598 //-----------------------------------------------------------------------------
599 
600 } // namespace basix
Definition: finite-element.h:191
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
cell::type cell_type() const
Definition: finite-element.cpp:415
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::vector< xt::xtensor< double, 2 >> &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)
Definition: finite-element.cpp:250
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
xt::xtensor< double, 4 > tabulate(int nd, const xt::xarray< double > &x) const
Definition: finite-element.cpp:457
xt::xtensor< double, 3 > map_pull_back(const xt::xtensor< double, 3 > &u, const xt::xtensor< double, 3 > &J, const tcb::span< const double > &detJ, const xt::xtensor< double, 3 > &K) const
Definition: finite-element.cpp:608
int value_size() const
Definition: finite-element.cpp:419
FiniteElement(FiniteElement &&element)=default
Move constructor.
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:436
void map_pull_back_m(const xt::xtensor< T, 3 > &u, const xt::xtensor< double, 3 > &J, const tcb::span< const double > &detJ, const xt::xtensor< double, 3 > &K, E &&U) const
Definition: finite-element.h:577
xt::xtensor< double, 3 > map_push_forward(const xt::xtensor< double, 3 > &U, const xt::xtensor< double, 3 > &J, const tcb::span< const double > &detJ, const xt::xtensor< double, 3 > &K) const
Definition: finite-element.cpp:596
const xt::xtensor< double, 2 > & points() const
Definition: finite-element.cpp:594
void unpermute_dofs(tcb::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:684
const xt::xtensor< double, 2 > & interpolation_matrix() const
Definition: finite-element.cpp:446
int degree() const
Definition: finite-element.cpp:417
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:425
~FiniteElement()=default
Destructor.
int num_points() const
Return the number of interpolation points.
Definition: finite-element.cpp:592
const std::vector< std::vector< int > > & entity_dofs() const
Definition: finite-element.cpp:451
int dim() const
Definition: finite-element.cpp:430
element::family family() const
Definition: finite-element.cpp:432
maps::type map_type
Element map type.
Definition: finite-element.h:462
void map_push_forward_m(const xt::xtensor< T, 3 > &U, const xt::xtensor< double, 3 > &J, const tcb::span< const double > &detJ, const xt::xtensor< double, 3 > &K, E &&u) const
Definition: finite-element.h:549
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:441
FiniteElement(const FiniteElement &element)=default
Copy constructor.
maps::type mapping_type() const
Definition: finite-element.cpp:434
void permute_dofs(tcb::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:618
xt::xtensor< double, 3 > base_transformations() const
Definition: finite-element.cpp:539
type
Cell type.
Definition: cell.h:21
void apply_map(O &&u, const P &U, const Mat0 &J, double detJ, const Mat1 &K, maps::type map_type)
Definition: maps.h:121
type
Cell type.
Definition: maps.h:24
Placeholder.
Definition: brezzi-douglas-marini.h:10
const char * cell_type(int handle)
Definition: basix.cpp:189
std::string version()
Definition: finite-element.cpp:750
FiniteElement create_element(std::string family, std::string cell, int degree)
Create an element by name.
Definition: finite-element.cpp:75
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)
Definition: finite-element.cpp:150
int degree(int handle)
Degree.
Definition: basix.cpp:195