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 "mappings.h"
13 #include "span.hpp"
14 #include <Eigen/Core>
15 #include <string>
16 #include <vector>
17 
19 namespace basix
20 {
21 
152 Eigen::MatrixXd compute_expansion_coefficients(
153  cell::type cell_type, const Eigen::MatrixXd& B, const Eigen::MatrixXd& M,
154  const Eigen::ArrayXXd& x, int degree, double kappa_tol = 0.0);
155 
171 std::pair<Eigen::ArrayXXd, Eigen::MatrixXd> combine_interpolation_data(
172  const Eigen::ArrayXXd& points_1d, const Eigen::ArrayXXd& points_2d,
173  const Eigen::ArrayXXd& points_3d, const Eigen::MatrixXd& matrix_1d,
174  const Eigen::MatrixXd& matrix_2d, const Eigen::MatrixXd& matrix_3d,
175  const int tdim, const int value_size);
176 
181 {
182 
183 public:
196  FiniteElement(element::family family, cell::type cell_type, int degree,
197  const std::vector<int>& value_shape,
198  const Eigen::ArrayXXd& coeffs,
199  const std::vector<std::vector<int>>& entity_dofs,
200  const std::vector<Eigen::MatrixXd>& base_transformations,
201  const Eigen::ArrayXXd& points, const Eigen::MatrixXd M = {},
202  mapping::type map_type = mapping::type::identity);
203 
205  FiniteElement(const FiniteElement& element) = default;
206 
208  FiniteElement(FiniteElement&& element) = default;
209 
211  ~FiniteElement() = default;
212 
214  FiniteElement& operator=(const FiniteElement& element) = default;
215 
217  FiniteElement& operator=(FiniteElement&& element) = default;
218 
233  std::vector<Eigen::ArrayXXd> tabulate(int nd, const Eigen::ArrayXXd& x) const;
234 
239  void tabulate(int nd, const Eigen::ArrayXXd& x, double* basis_data) const;
240 
243  cell::type cell_type() const;
244 
247  int degree() const;
248 
252  int value_size() const;
253 
256  const std::vector<int>& value_shape() const;
257 
261  int dim() const;
262 
265  element::family family() const;
266 
269  const mapping::type mapping_type() const;
270 
277  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
278  map_push_forward(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
279  Eigen::RowMajor>& reference_data,
280  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
281  Eigen::RowMajor>& J,
282  const tcb::span<const double>& detJ,
283  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
284  Eigen::RowMajor>& K) const;
285 
292  template <typename T>
293  void
294  map_push_forward_m(const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic,
295  Eigen::RowMajor>& reference_data,
296  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
297  Eigen::RowMajor>& J,
298  const tcb::span<const double>& detJ,
299  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
300  Eigen::RowMajor>& K,
301  T* physical_data) const;
302 
309  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
310  map_pull_back(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
311  Eigen::RowMajor>& physical_data,
312  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
313  Eigen::RowMajor>& J,
314  const tcb::span<const double>& detJ,
315  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
316  Eigen::RowMajor>& K) const;
317 
324  template <typename T>
325  void map_pull_back_m(
326  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>& physical_data,
327  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
328  Eigen::RowMajor>& J,
329  const tcb::span<const double>& detJ,
330  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic,
331  Eigen::RowMajor>& K,
332  T* reference_data) const;
333 
340  const std::vector<std::vector<int>>& entity_dofs() const;
341 
420  std::vector<Eigen::MatrixXd> base_permutations() const;
421 
423  const Eigen::ArrayXXd& points() const;
424 
426  int num_points() const;
427 
434  const Eigen::MatrixXd& interpolation_matrix() const;
435 
436 private:
437  static int compute_value_size(mapping::type map_type, int dim);
438 
439  // Cell type
440  cell::type _cell_type;
441 
442  // The name of the finite element family
443  element::family _family;
444 
445  // Degree
446  int _degree;
447 
448  // Value shape
449  std::vector<int> _value_shape;
450 
452  mapping::type _map_type;
453 
454  // Shape function coefficient of expansion sets on cell. If shape
455  // function is given by @f$\psi_i = \sum_{k} \phi_{k}
456  // \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. i.e.,
457  // _coeffs.row(i) are the expansion coefficients for shape function i
458  // (@f$\psi_{i}@f$).
459  Eigen::MatrixXd _coeffs;
460 
461  // Number of dofs associated each subentity
462  //
463  // The dofs of an element are associated with entities of different
464  // topological dimension (vertices, edges, faces, cells). The dofs are
465  // listed in this order, with vertex dofs first. Each entry is the dof
466  // count on the associated entity, as listed by cell::topology.
467  std::vector<std::vector<int>> _entity_dofs;
468 
469  // Base transformations
470  std::vector<Eigen::MatrixXd> _base_transformations;
471 
472  // Set of points used for point evaluation
473  // Experimental - currently used for an implementation of
474  // "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
475  // away. For non-Lagrange elements, these points will be used in combination
476  // with _interpolation_matrix to perform interpolation
477  Eigen::ArrayXXd _points;
478 
480  Eigen::MatrixXd _matM;
481 
482  // The mapping that maps values on the reference to values on a physical cell
483  std::function<std::vector<double>(const tcb::span<const double>&,
484  const Eigen::MatrixXd&, const double,
485  const Eigen::MatrixXd&)>
486  _map_push_forward;
487 };
488 
489 template <typename T>
491  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
492  reference_data,
493  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
494  J,
495  const tcb::span<const double>& detJ,
496  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
497  K,
498  T* physical_data) const
499 {
500  const int reference_dim = cell::topological_dimension(_cell_type);
501  const int physical_dim = J.cols() / reference_dim;
502  const int physical_value_size = compute_value_size(_map_type, physical_dim);
503  const int reference_value_size = value_size();
504  const int npoints = J.rows();
505  const int nresults = reference_data.rows() / npoints;
506 
507  for (int pt = 0; pt < npoints; ++pt)
508  {
509  Eigen::Map<
510  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>>
511  reference_block(reference_data.row(pt).data(), reference_value_size,
512  nresults);
513  Eigen::Map<Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>>
514  physical_block(physical_data + pt * physical_value_size * nresults,
515  physical_value_size, nresults);
516  Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
517  Eigen::RowMajor>>
518  current_J(J.row(pt).data(), physical_dim, reference_dim);
519  Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
520  Eigen::RowMajor>>
521  current_K(K.row(pt).data(), reference_dim, physical_dim);
522  if constexpr (std::is_same<T, double>::value)
523  {
524  for (int i = 0; i < reference_block.cols(); ++i)
525  {
526  Eigen::ArrayXd col = reference_block.col(i);
527  std::vector<double> u
528  = _map_push_forward(col, current_J, detJ[pt], current_K);
529  for (std::size_t j = 0; j < u.size(); ++j)
530  physical_block(j, i) = u[j];
531  }
532  }
533  else
534  {
535  for (int i = 0; i < reference_block.cols(); ++i)
536  {
537  Eigen::ArrayXd tmp_r = reference_block.col(i).real();
538  Eigen::ArrayXd tmp_c = reference_block.col(i).imag();
539  std::vector<double> ur
540  = _map_push_forward(tmp_r, current_J, detJ[pt], current_K);
541  std::vector<double> uc
542  = _map_push_forward(tmp_c, current_J, detJ[pt], current_K);
543  for (std::size_t j = 0; j < ur.size(); ++j)
544  physical_block(j, i) = std::complex(ur[j], uc[j]);
545  }
546  }
547  }
548 }
549 //-----------------------------------------------------------------------------
550 template <typename T>
552  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>& physical_data,
553  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
554  J,
555  const tcb::span<const double>& detJ,
556  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
557  K,
558  T* reference_data) const
559 {
560  const int reference_dim = cell::topological_dimension(_cell_type);
561  const int physical_dim = J.cols() / reference_dim;
562  const int reference_value_size = value_size();
563  const int npoints = J.rows();
564  const int nresults = physical_data.rows() / npoints;
565 
566  Eigen::Map<Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>>
567  reference_array(reference_data, nresults * npoints, reference_value_size);
568 
569  for (int pt = 0; pt < npoints; ++pt)
570  {
571  Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
572  Eigen::RowMajor>>
573  current_J(J.row(pt).data(), physical_dim, reference_dim);
574  Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
575  Eigen::RowMajor>>
576  current_K(K.row(pt).data(), reference_dim, physical_dim);
577  if constexpr (std::is_same<T, double>::value)
578  {
579  for (int i = 0; i < nresults; ++i)
580  {
581  for (int i = 0; i < nresults; ++i)
582  {
583  Eigen::ArrayXd tmp = physical_data.row(pt * nresults + i);
584  std::vector<double> U
585  = _map_push_forward(tmp, current_K, 1 / detJ[pt], current_J);
586  for (std::size_t j = 0; j < U.size(); ++j)
587  reference_array(pt * nresults + i, j) = U[j];
588  }
589  }
590  }
591  else
592  {
593  for (int i = 0; i < nresults; ++i)
594  {
595  for (int i = 0; i < nresults; ++i)
596  {
597  Eigen::ArrayXd tmp_r = physical_data.row(pt * nresults + i).real();
598  Eigen::ArrayXd tmp_c = physical_data.row(pt * nresults + i).imag();
599  std::vector<double> Ur
600  = _map_push_forward(tmp_r, current_K, 1 / detJ[pt], current_J);
601  std::vector<double> Uc
602  = _map_push_forward(tmp_c, current_K, 1 / detJ[pt], current_J);
603  for (std::size_t j = 0; j < Ur.size(); ++j)
604  reference_array(pt * nresults + i, j) = std::complex(Ur[j], Uc[j]);
605  }
606  }
607  }
608  }
609 }
610 //-----------------------------------------------------------------------------
611 
613 FiniteElement create_element(std::string family, std::string cell, int degree);
614 
616 FiniteElement create_element(element::family family, cell::type cell,
617  int degree);
618 
621 std::string version();
622 
623 } // namespace basix
Definition: finite-element.h:181
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
cell::type cell_type() const
Definition: finite-element.cpp:179
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
const Eigen::ArrayXXd & points() const
Return a set of interpolation points.
Definition: finite-element.cpp:271
Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > map_push_forward(const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &reference_data, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &J, const tcb::span< const double > &detJ, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &K) const
Definition: finite-element.cpp:274
void map_push_forward_m(const Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &reference_data, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &J, const tcb::span< const double > &detJ, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &K, T *physical_data) const
Definition: finite-element.h:490
const mapping::type mapping_type() const
Definition: finite-element.cpp:200
int value_size() const
Definition: finite-element.cpp:183
FiniteElement(FiniteElement &&element)=default
Move constructor.
void map_pull_back_m(const Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic > &physical_data, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &J, const tcb::span< const double > &detJ, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &K, T *reference_data) const
Definition: finite-element.h:551
Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > map_pull_back(const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &physical_data, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &J, const tcb::span< const double > &detJ, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > &K) const
Definition: finite-element.cpp:322
int degree() const
Definition: finite-element.cpp:181
const std::vector< int > & value_shape() const
Definition: finite-element.cpp:191
~FiniteElement()=default
Destructor.
int num_points() const
Return the number of interpolation points.
Definition: finite-element.cpp:269
const std::vector< std::vector< int > > & entity_dofs() const
Definition: finite-element.cpp:207
int dim() const
Definition: finite-element.cpp:196
element::family family() const
Definition: finite-element.cpp:198
const Eigen::MatrixXd & interpolation_matrix() const
Definition: finite-element.cpp:202
FiniteElement(const FiniteElement &element)=default
Copy constructor.
std::vector< Eigen::ArrayXXd > tabulate(int nd, const Eigen::ArrayXXd &x) const
Definition: finite-element.cpp:213
FiniteElement(element::family family, cell::type cell_type, int degree, const std::vector< int > &value_shape, const Eigen::ArrayXXd &coeffs, const std::vector< std::vector< int >> &entity_dofs, const std::vector< Eigen::MatrixXd > &base_transformations, const Eigen::ArrayXXd &points, const Eigen::MatrixXd M={}, mapping::type map_type=mapping::type::identity)
Definition: finite-element.cpp:154
std::vector< Eigen::MatrixXd > base_permutations() const
Definition: finite-element.cpp:264
type
Cell type.
Definition: cell.h:21
int topological_dimension(cell::type celltype)
Definition: cell.cpp:136
type
Cell type.
Definition: mappings.h:18
Placeholder.
Definition: basix.h:10
Eigen::MatrixXd compute_expansion_coefficients(cell::type cell_type, const Eigen::MatrixXd &B, const Eigen::MatrixXd &M, const Eigen::ArrayXXd &x, int degree, double kappa_tol=0.0)
Definition: finite-element.cpp:83
const char * cell_type(int handle)
Definition: basix.cpp:161
std::string version()
Definition: finite-element.cpp:369
FiniteElement create_element(std::string family, std::string cell, int degree)
Create an element by name.
Definition: finite-element.cpp:24
std::pair< Eigen::ArrayXXd, Eigen::MatrixXd > combine_interpolation_data(const Eigen::ArrayXXd &points_1d, const Eigen::ArrayXXd &points_2d, const Eigen::ArrayXXd &points_3d, const Eigen::MatrixXd &matrix_1d, const Eigen::MatrixXd &matrix_2d, const Eigen::MatrixXd &matrix_3d, const int tdim, const int value_size)
Definition: finite-element.cpp:118
int degree(int handle)
Degree.
Definition: basix.cpp:167