Home Installation C++ docs Python docs
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 "lattice.h"
13 #include "maps.h"
14 #include "precompute.h"
15 #include <array>
16 #include <numeric>
17 #include <set>
18 #include <string>
19 #include <vector>
20 #include <xtensor/xadapt.hpp>
21 #include <xtensor/xcomplex.hpp>
22 #include <xtensor/xtensor.hpp>
23 #include <xtensor/xview.hpp>
24 #include <xtl/xspan.hpp>
25 
27 namespace basix
28 {
29 
185 xt::xtensor<double, 3> compute_expansion_coefficients(
186  cell::type cell_type, const xt::xtensor<double, 2>& B,
187  const std::vector<std::vector<xt::xtensor<double, 3>>>& M,
188  const std::vector<std::vector<xt::xtensor<double, 2>>>& x, int degree,
189  double kappa_tol = 0.0);
190 
195 {
196 
197 public:
212  FiniteElement(element::family family, cell::type cell_type, int degree,
213  const std::vector<std::size_t>& value_shape,
214  const xt::xtensor<double, 3>& coeffs,
215  const std::map<cell::type, xt::xtensor<double, 3>>&
217  const std::array<std::vector<xt::xtensor<double, 2>>, 4>& x,
218  const std::array<std::vector<xt::xtensor<double, 3>>, 4>& M,
219  maps::type map_type = maps::type::identity);
220 
222  FiniteElement(const FiniteElement& element) = default;
223 
225  FiniteElement(FiniteElement&& element) = default;
226 
228  ~FiniteElement() = default;
229 
231  FiniteElement& operator=(const FiniteElement& element) = default;
232 
234  FiniteElement& operator=(FiniteElement&& element) = default;
235 
253  xt::xtensor<double, 4> tabulate(int nd, const xt::xarray<double>& x) const;
254 
259  void tabulate(int nd, const xt::xarray<double>& x,
260  xt::xtensor<double, 4>& basis_data) const;
261 
264  cell::type cell_type() const;
265 
268  int degree() const;
269 
273  int value_size() const;
274 
277  const std::vector<int>& value_shape() const;
278 
282  int dim() const;
283 
286  element::family family() const;
287 
290  maps::type mapping_type() const;
291 
295 
299 
314  xt::xtensor<double, 3>
315  map_push_forward(const xt::xtensor<double, 3>& U,
316  const xt::xtensor<double, 3>& J,
317  const xtl::span<const double>& detJ,
318  const xt::xtensor<double, 3>& K) const;
319 
335  template <typename O, typename P, typename Q, typename S, typename T>
336  void map_push_forward_m(const O& U, const P& J, const Q& detJ, const S& K,
337  T&& u) const
338  {
339  // FIXME: Should U.shape(2) be replaced by the physical value size?
340  // Can it differ?
341 
342  // Loop over points that share J
343  for (std::size_t i = 0; i < U.shape(0); ++i)
344  {
345  auto _J = xt::view(J, i, xt::all(), xt::all());
346  auto _K = xt::view(K, i, xt::all(), xt::all());
347  auto _U = xt::view(U, i, xt::all(), xt::all());
348  auto _u = xt::view(u, i, xt::all(), xt::all());
349  maps::apply_map(_u, _U, _J, detJ[i], _K, map_type);
350  }
351  }
352 
359  xt::xtensor<double, 3> map_pull_back(const xt::xtensor<double, 3>& u,
360  const xt::xtensor<double, 3>& J,
361  const xtl::span<const double>& detJ,
362  const xt::xtensor<double, 3>& K) const;
363 
379  template <typename O, typename P, typename Q, typename S, typename T>
380  void map_pull_back_m(const O& u, const P& J, const Q& detJ, const S& K,
381  T&& U) const
382  {
383  // Loop over points that share K and K
384  for (std::size_t i = 0; i < u.shape(0); ++i)
385  {
386  auto _J = xt::view(J, i, xt::all(), xt::all());
387  auto _K = xt::view(K, i, xt::all(), xt::all());
388  auto _u = xt::view(u, i, xt::all(), xt::all());
389  auto _U = xt::view(U, i, xt::all(), xt::all());
390  maps::apply_map(_U, _u, _K, 1.0 / detJ[i], _J, map_type);
391  }
392  }
393 
406  const std::vector<std::vector<int>>& num_entity_dofs() const;
407 
413  const std::vector<std::vector<int>>& num_entity_closure_dofs() const;
414 
421  const std::vector<std::vector<std::set<int>>>& entity_dofs() const;
422 
429  const std::vector<std::vector<std::set<int>>>& entity_closure_dofs() const;
430 
509  xt::xtensor<double, 3> base_transformations() const;
510 
512  std::map<cell::type, xt::xtensor<double, 3>> entity_transformations() const;
513 
517  void permute_dofs(const xtl::span<std::int32_t>& dofs,
518  std::uint32_t cell_info) const;
519 
523  void unpermute_dofs(const xtl::span<std::int32_t>& dofs,
524  std::uint32_t cell_info) const;
525 
530  template <typename T>
531  void apply_dof_transformation(const xtl::span<T>& data, int block_size,
532  std::uint32_t cell_info) const;
533 
538  template <typename T>
539  void apply_transpose_dof_transformation(const xtl::span<T>& data,
540  int block_size,
541  std::uint32_t cell_info) const;
542 
547  template <typename T>
549  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
550 
555  template <typename T>
556  void apply_inverse_dof_transformation(const xtl::span<T>& data,
557  int block_size,
558  std::uint32_t cell_info) const;
559 
564  template <typename T>
565  void apply_dof_transformation_to_transpose(const xtl::span<T>& data,
566  int block_size,
567  std::uint32_t cell_info) const;
568 
573  template <typename T>
575  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
576 
581  template <typename T>
583  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
584 
589  template <typename T>
591  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const;
592 
597  const xt::xtensor<double, 2>& points() const;
598 
600  int num_points() const;
601 
608  const xt::xtensor<double, 2>& interpolation_matrix() const;
609 
616  template <typename T>
617  void interpolate(const xtl::span<T>& coefficients,
618  const xtl::span<const T>& data, const int block_size) const;
619 
622 
623 private:
624  // Cell type
625  cell::type _cell_type;
626 
627  // Topological dimension of the cell
628  std::size_t _cell_tdim;
629 
630  // Topological dimension of the cell
631  std::vector<std::vector<cell::type>> _cell_subentity_types;
632 
633  // Finite element family
634  element::family _family;
635 
636  // Degree
637  int _degree;
638 
639  // Value shape
640  std::vector<int> _value_shape;
641 
643  maps::type _map_type;
644 
645  // Shape function coefficient of expansion sets on cell. If shape
646  // function is given by @f$\psi_i = \sum_{k} \phi_{k}
647  // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. i.e.,
648  // _coeffs.row(i) are the expansion coefficients for shape function i
649  // (@f$\psi_{i}@f$).
650  xt::xtensor<double, 2> _coeffs;
651 
652  // Number of dofs associated with each cell (sub-)entity
653  //
654  // The dofs of an element are associated with entities of different
655  // topological dimension (vertices, edges, faces, cells). The dofs are
656  // listed in this order, with vertex dofs first. Each entry is the dof
657  // count on the associated entity, as listed by cell::topology.
658  std::vector<std::vector<int>> _num_edofs;
659 
660  // Number of dofs associated with the closure of each cell (sub-)entity
661  std::vector<std::vector<int>> _num_e_closure_dofs;
662 
663  // Dofs associated with each cell (sub-)entity
664  std::vector<std::vector<std::set<int>>> _edofs;
665 
666  // Dofs associated with each cell (sub-)entity
667  std::vector<std::vector<std::set<int>>> _e_closure_dofs;
668 
669  // Entity transformations
670  std::map<cell::type, xt::xtensor<double, 3>> _entity_transformations;
671 
672  // Set of points used for point evaluation
673  // Experimental - currently used for an implementation of
674  // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
675  // away. For non-Lagrange elements, these points will be used in combination
676  // with _interpolation_matrix to perform interpolation
677  xt::xtensor<double, 2> _points;
678 
679  // Interpolation points on the cell. The shape is (entity_dim, num
680  // entities of given dimension, num_points, tdim)
681  std::array<std::vector<xt::xtensor<double, 2>>, 4> _x;
682 
684  xt::xtensor<double, 2> _matM;
685 
687  std::array<std::vector<xt::xtensor<double, 3>>, 4> _matM_new;
688 
690  bool _dof_transformations_are_permutations;
691 
693  bool _dof_transformations_are_identity;
694 
698  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
699 
703  std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_rev;
704 
706  std::map<cell::type,
707  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
708  xt::xtensor<double, 2>>>>
709  _etrans;
710 
712  std::map<cell::type,
713  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
714  xt::xtensor<double, 2>>>>
715  _etransT;
716 
718  std::map<cell::type,
719  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
720  xt::xtensor<double, 2>>>>
721  _etrans_inv;
722 
724  std::map<cell::type,
725  std::vector<std::tuple<std::vector<std::size_t>, std::vector<double>,
726  xt::xtensor<double, 2>>>>
727  _etrans_invT;
728 };
729 
737 create_element(element::family family, cell::type cell, int degree,
738  lattice::type lattice_type);
739 
745 create_element(element::family family, cell::type cell, int degree);
746 
749 std::string version();
750 
751 //-----------------------------------------------------------------------------
752 template <typename T>
753 void FiniteElement::apply_dof_transformation(const xtl::span<T>& data,
754  int block_size,
755  std::uint32_t cell_info) const
756 {
757  if (_dof_transformations_are_identity)
758  return;
759 
760  if (_cell_tdim >= 2)
761  {
762  // This assumes 3 bits are used per face. This will need updating if
763  // 3D cells with faces with more than 4 sides are implemented
764  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
765  int dofstart
766  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
767 
768  // Transform DOFs on edges
769  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
770  {
771  // Reverse an edge
772  if (cell_info >> (face_start + e) & 1)
773  precompute::apply_matrix(_etrans.at(cell::type::interval)[0], data,
774  dofstart, block_size);
775  dofstart += _num_edofs[1][e];
776  }
777 
778  if (_cell_tdim == 3)
779  {
780  // Permute DOFs on faces
781  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
782  {
783  // Reflect a face
784  if (cell_info >> (3 * f) & 1)
785  precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[1],
786  data, dofstart, block_size);
787 
788  // Rotate a face
789  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
790  precompute::apply_matrix(_etrans.at(_cell_subentity_types[2][f])[0],
791  data, dofstart, block_size);
792  dofstart += _num_edofs[2][f];
793  }
794  }
795  }
796 }
797 //-----------------------------------------------------------------------------
798 template <typename T>
800  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
801 {
802  if (_dof_transformations_are_identity)
803  return;
804 
805  if (_cell_tdim >= 2)
806  {
807  // This assumes 3 bits are used per face. This will need updating if
808  // 3D cells with faces with more than 4 sides are implemented
809  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
810  int dofstart
811  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
812 
813  // Transform DOFs on edges
814  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
815  {
816  // Reverse an edge
817  if (cell_info >> (face_start + e) & 1)
818  precompute::apply_matrix(_etransT.at(cell::type::interval)[0], data,
819  dofstart, block_size);
820  dofstart += _num_edofs[1][e];
821  }
822 
823  if (_cell_tdim == 3)
824  {
825  // Permute DOFs on faces
826  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
827  {
828  // Rotate a face
829  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
830  precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[0],
831  data, dofstart, block_size);
832  // Reflect a face
833  if (cell_info >> (3 * f) & 1)
834  precompute::apply_matrix(_etransT.at(_cell_subentity_types[2][f])[1],
835  data, dofstart, block_size);
836  dofstart += _num_edofs[2][f];
837  }
838  }
839  }
840 }
841 //-----------------------------------------------------------------------------
842 template <typename T>
844  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
845 {
846  if (_dof_transformations_are_identity)
847  return;
848 
849  if (_cell_tdim >= 2)
850  {
851  // This assumes 3 bits are used per face. This will need updating if 3D
852  // cells with faces with more than 4 sides are implemented
853  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
854  int dofstart
855  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
856 
857  // Transform DOFs on edges
858  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
859  {
860  // Reverse an edge
861  if (cell_info >> (face_start + e) & 1)
862  precompute::apply_matrix(_etrans_invT.at(cell::type::interval)[0], data,
863  dofstart, block_size);
864  dofstart += _num_edofs[1][e];
865  }
866 
867  if (_cell_tdim == 3)
868  {
869  // Permute DOFs on faces
870  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
871  {
872  // Reflect a face
873  if (cell_info >> (3 * f) & 1)
875  _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
876  block_size);
877 
878  // Rotate a face
879  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
881  _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
882  block_size);
883  dofstart += _num_edofs[2][f];
884  }
885  }
886  }
887 }
888 //-----------------------------------------------------------------------------
889 template <typename T>
891  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
892 {
893  if (_dof_transformations_are_identity)
894  return;
895 
896  if (_cell_tdim >= 2)
897  {
898  // This assumes 3 bits are used per face. This will need updating if 3D
899  // cells with faces with more than 4 sides are implemented
900  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
901  int dofstart
902  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
903 
904  // Transform DOFs on edges
905  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
906  {
907  // Reverse an edge
908  if (cell_info >> (face_start + e) & 1)
909  precompute::apply_matrix(_etrans_inv.at(cell::type::interval)[0], data,
910  dofstart, block_size);
911  dofstart += _num_edofs[1][e];
912  }
913 
914  if (_cell_tdim == 3)
915  {
916  // Permute DOFs on faces
917  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
918  {
919  // Rotate a face
920  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
922  _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
923  block_size);
924  // Reflect a face
925  if (cell_info >> (3 * f) & 1)
927  _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
928  block_size);
929  dofstart += _num_edofs[2][f];
930  }
931  }
932  }
933 }
934 //-----------------------------------------------------------------------------
935 template <typename T>
937  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
938 {
939  if (_dof_transformations_are_identity)
940  return;
941 
942  if (_cell_tdim >= 2)
943  {
944  // This assumes 3 bits are used per face. This will need updating if
945  // 3D cells with faces with more than 4 sides are implemented
946  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
947  int dofstart
948  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
949 
950  // Transform DOFs on edges
951  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
952  {
953  // Reverse an edge
954  if (cell_info >> (face_start + e) & 1)
956  _etrans.at(cell::type::interval)[0], data, dofstart, block_size);
957  dofstart += _num_edofs[1][e];
958  }
959 
960  if (_cell_tdim == 3)
961  {
962  // Permute DOFs on faces
963  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
964  {
965  // Reflect a face
966  if (cell_info >> (3 * f) & 1)
968  _etrans.at(_cell_subentity_types[2][f])[1], data, dofstart,
969  block_size);
970 
971  // Rotate a face
972  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
974  _etrans.at(_cell_subentity_types[2][f])[0], data, dofstart,
975  block_size);
976  dofstart += _num_edofs[2][f];
977  }
978  }
979  }
980 }
981 //-----------------------------------------------------------------------------
982 template <typename T>
984  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
985 {
986  if (_dof_transformations_are_identity)
987  return;
988 
989  if (_cell_tdim >= 2)
990  {
991  // This assumes 3 bits are used per face. This will need updating if
992  // 3D cells with faces with more than 4 sides are implemented
993  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
994  int dofstart
995  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
996 
997  // Transform DOFs on edges
998  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
999  {
1000  // Reverse an edge
1001  if (cell_info >> (face_start + e) & 1)
1003  _etrans_invT.at(cell::type::interval)[0], data, dofstart,
1004  block_size);
1005  dofstart += _num_edofs[1][e];
1006  }
1007 
1008  if (_cell_tdim == 3)
1009  {
1010  // Permute DOFs on faces
1011  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1012  {
1013  // Reflect a face
1014  if (cell_info >> (3 * f) & 1)
1016  _etrans_invT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1017  block_size);
1018 
1019  // Rotate a face
1020  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1022  _etrans_invT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1023  block_size);
1024  dofstart += _num_edofs[2][f];
1025  }
1026  }
1027  }
1028 }
1029 //-----------------------------------------------------------------------------
1030 template <typename T>
1032  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1033 {
1034  if (_dof_transformations_are_identity)
1035  return;
1036 
1037  if (_cell_tdim >= 2)
1038  {
1039  // This assumes 3 bits are used per face. This will need updating if
1040  // 3D cells with faces with more than 4 sides are implemented
1041  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1042  int dofstart
1043  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1044 
1045  // Transform DOFs on edges
1046  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1047  {
1048  // Reverse an edge
1049  if (cell_info >> (face_start + e) & 1)
1051  _etransT.at(cell::type::interval)[0], data, dofstart, block_size);
1052  dofstart += _num_edofs[1][e];
1053  }
1054 
1055  if (_cell_tdim == 3)
1056  {
1057  // Permute DOFs on faces
1058  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1059  {
1060  // Rotate a face
1061  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1063  _etransT.at(_cell_subentity_types[2][f])[0], data, dofstart,
1064  block_size);
1065 
1066  // Reflect a face
1067  if (cell_info >> (3 * f) & 1)
1069  _etransT.at(_cell_subentity_types[2][f])[1], data, dofstart,
1070  block_size);
1071 
1072  dofstart += _num_edofs[2][f];
1073  }
1074  }
1075  }
1076 }
1077 //-----------------------------------------------------------------------------
1078 template <typename T>
1080  const xtl::span<T>& data, int block_size, std::uint32_t cell_info) const
1081 {
1082  if (_dof_transformations_are_identity)
1083  return;
1084 
1085  if (_cell_tdim >= 2)
1086  {
1087  // This assumes 3 bits are used per face. This will need updating if
1088  // 3D cells with faces with more than 4 sides are implemented
1089  int face_start = _cell_tdim == 3 ? 3 * _num_edofs[2].size() : 0;
1090  int dofstart
1091  = std::accumulate(_num_edofs[0].cbegin(), _num_edofs[0].cend(), 0);
1092 
1093  // Transform DOFs on edges
1094  for (std::size_t e = 0; e < _num_edofs[1].size(); ++e)
1095  {
1096  // Reverse an edge
1097  if (cell_info >> (face_start + e) & 1)
1099  _etrans_inv.at(cell::type::interval)[0], data, dofstart,
1100  block_size);
1101  dofstart += _num_edofs[1][e];
1102  }
1103 
1104  if (_cell_tdim == 3)
1105  {
1106  // Permute DOFs on faces
1107  for (std::size_t f = 0; f < _num_edofs[2].size(); ++f)
1108  {
1109  // Rotate a face
1110  for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
1112  _etrans_inv.at(_cell_subentity_types[2][f])[0], data, dofstart,
1113  block_size);
1114 
1115  // Reflect a face
1116  if (cell_info >> (3 * f) & 1)
1118  _etrans_inv.at(_cell_subentity_types[2][f])[1], data, dofstart,
1119  block_size);
1120 
1121  dofstart += _num_edofs[2][f];
1122  }
1123  }
1124  }
1125 }
1126 //-----------------------------------------------------------------------------
1127 template <typename T>
1128 void FiniteElement::interpolate(const xtl::span<T>& coefficients,
1129  const xtl::span<const T>& data,
1130  const int block_size) const
1131 {
1132  if (block_size != 1)
1133  {
1134  throw std::runtime_error(
1135  "Interpolation of blocked data not implemented yet.");
1136  }
1137 
1138  const std::size_t rows = dim();
1139 
1140  // Compute coefficients = Pi * x (matrix-vector multiply)
1141  const xt::xtensor<double, 2>& Pi = interpolation_matrix();
1142  assert(Pi.size() % rows == 0);
1143  const std::size_t cols = Pi.size() / rows;
1144  for (std::size_t i = 0; i < rows; ++i)
1145  {
1146  // Can be replaced with std::transform_reduce once GCC 8 series dies.
1147  // Dot product between row i of the matrix and 'data'
1148  coefficients[i] = std::inner_product(std::next(Pi.data(), i * cols),
1149  std::next(Pi.data(), i * cols + cols),
1150  data.data(), T(0.0));
1151  }
1152 }
1153 //-----------------------------------------------------------------------------
1154 
1155 } // namespace basix
basix::FiniteElement::permute_dofs
void permute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:678
basix::FiniteElement::entity_dofs
const std::vector< std::vector< std::set< int > > > & entity_dofs() const
Definition: finite-element.cpp:530
basix::FiniteElement::dof_transformations_are_identity
bool dof_transformations_are_identity() const
Definition: finite-element.cpp:514
basix::FiniteElement::tabulate
xt::xtensor< double, 4 > tabulate(int nd, const xt::xarray< double > &x) const
Definition: finite-element.cpp:548
basix::FiniteElement::apply_transpose_dof_transformation
void apply_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:799
basix::FiniteElement::apply_dof_transformation
void apply_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:753
basix::FiniteElement::degree
int degree() const
Definition: finite-element.cpp:490
basix::FiniteElement::~FiniteElement
~FiniteElement()=default
Destructor.
basix::FiniteElement::num_entity_dofs
const std::vector< std::vector< int > > & num_entity_dofs() const
Definition: finite-element.cpp:524
basix::cell::type
type
Cell type.
Definition: cell.h:16
basix::FiniteElement::points
const xt::xtensor< double, 2 > & points() const
Definition: finite-element.cpp:655
basix::FiniteElement::entity_transformations
std::map< cell::type, xt::xtensor< double, 3 > > entity_transformations() const
Return the entity dof transformation matricess.
Definition: finite-element.cpp:780
basix::FiniteElement::cell_type
cell::type cell_type() const
Definition: finite-element.cpp:488
basix::FiniteElement::apply_inverse_transpose_dof_transformation_to_transpose
void apply_inverse_transpose_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:983
basix::FiniteElement::apply_inverse_dof_transformation
void apply_inverse_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:890
basix::FiniteElement::mapping_type
maps::type mapping_type() const
Definition: finite-element.cpp:507
basix::FiniteElement::apply_transpose_dof_transformation_to_transpose
void apply_transpose_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1031
basix::FiniteElement::dof_transformations_are_permutations
bool dof_transformations_are_permutations() const
Definition: finite-element.cpp:509
basix::compute_expansion_coefficients
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:159
basix::FiniteElement::family
element::family family() const
Definition: finite-element.cpp:505
basix::FiniteElement::apply_inverse_dof_transformation_to_transpose
void apply_inverse_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:1079
basix::FiniteElement::dim
int dim() const
Definition: finite-element.cpp:503
basix::FiniteElement::map_pull_back
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
Definition: finite-element.cpp:668
basix::FiniteElement::apply_dof_transformation_to_transpose
void apply_dof_transformation_to_transpose(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:936
basix::maps::apply_map
void apply_map(O &&u, const P &U, const Mat0 &J, double detJ, const Mat1 &K, maps::type map_type)
Definition: maps.h:130
basix::create_element
FiniteElement create_element(element::family family, cell::type cell, int degree, lattice::type lattice_type)
Definition: finite-element.cpp:146
basix::maps::type
type
Cell type.
Definition: maps.h:19
basix::FiniteElement::value_shape
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:498
basix::FiniteElement::entity_closure_dofs
const std::vector< std::vector< std::set< int > > > & entity_closure_dofs() const
Definition: finite-element.cpp:542
basix::FiniteElement::base_transformations
xt::xtensor< double, 3 > base_transformations() const
Definition: finite-element.cpp:595
basix::FiniteElement::num_entity_closure_dofs
const std::vector< std::vector< int > > & num_entity_closure_dofs() const
Definition: finite-element.cpp:536
basix::precompute::apply_matrix_to_transpose
void apply_matrix_to_transpose(const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:340
basix::FiniteElement
Definition: finite-element.h:194
basix::FiniteElement::interpolation_matrix
const xt::xtensor< double, 2 > & interpolation_matrix() const
Definition: finite-element.cpp:519
basix::FiniteElement::map_push_forward
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
Definition: finite-element.cpp:657
basix::version
std::string version()
Definition: finite-element.cpp:785
basix::FiniteElement::value_size
int value_size() const
Definition: finite-element.cpp:492
basix::FiniteElement::map_type
maps::type map_type
Element map type.
Definition: finite-element.h:621
basix::FiniteElement::interpolate
void interpolate(const xtl::span< T > &coefficients, const xtl::span< const T > &data, const int block_size) const
Definition: finite-element.h:1128
basix::FiniteElement::map_push_forward_m
void map_push_forward_m(const O &U, const P &J, const Q &detJ, const S &K, T &&u) const
Definition: finite-element.h:336
basix
Placeholder.
Definition: brezzi-douglas-marini.h:9
basix::FiniteElement::unpermute_dofs
void unpermute_dofs(const xtl::span< std::int32_t > &dofs, std::uint32_t cell_info) const
Definition: finite-element.cpp:729
basix::precompute::apply_matrix
void apply_matrix(const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
Definition: precompute.h:311
basix::FiniteElement::num_points
int num_points() const
Return the number of interpolation points.
Definition: finite-element.cpp:653
basix::FiniteElement::operator=
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
basix::FiniteElement::map_pull_back_m
void map_pull_back_m(const O &u, const P &J, const Q &detJ, const S &K, T &&U) const
Definition: finite-element.h:380
basix::FiniteElement::apply_inverse_transpose_dof_transformation
void apply_inverse_transpose_dof_transformation(const xtl::span< T > &data, int block_size, std::uint32_t cell_info) const
Definition: finite-element.h:843
basix::FiniteElement::FiniteElement
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)
Definition: finite-element.cpp:259