Loading...
Searching...
No Matches
finite-element.h
// Copyright (c) 2020-2024 Chris Richardson, Matthew Scroggs and Garth . Wells
// FEniCS Project
// SPDX-License-Identifier: MIT
#pragma once
#include "cell.h"
#include "element-families.h"
#include "maps.h"
#include "mdspan.hpp"
#include "polyset.h"
#include "precompute.h"
#include "sobolev-spaces.h"
#include <array>
#include <concepts>
#include <cstdint>
#include <functional>
#include <map>
#include <numeric>
#include <span>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
namespace basix
{
namespace impl
{
template <typename T, std::size_t d>
using mdspan_t = md::mdspan<T, md::dextents<std::size_t, d>>;
template <typename T, std::size_t d>
using mdarray_t
= md::MDSPAN_IMPL_PROPOSED_NAMESPACE::mdarray<T,
md::dextents<std::size_t, d>>;
template <typename T>
std::array<std::vector<mdspan_t<const T, 2>>, 4>
to_mdspan(std::array<std::vector<mdarray_t<T, 2>>, 4>& x)
{
std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
for (std::size_t i = 0; i < x.size(); ++i)
for (std::size_t j = 0; j < x[i].size(); ++j)
x1[i].emplace_back(x[i][j].data(), x[i][j].extents());
return x1;
}
template <typename T>
std::array<std::vector<mdspan_t<const T, 4>>, 4>
to_mdspan(std::array<std::vector<mdarray_t<T, 4>>, 4>& M)
{
std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
for (std::size_t i = 0; i < M.size(); ++i)
for (std::size_t j = 0; j < M[i].size(); ++j)
M1[i].emplace_back(M[i][j].data(), M[i][j].extents());
return M1;
}
template <typename T>
std::array<std::vector<mdspan_t<const T, 2>>, 4>
to_mdspan(const std::array<std::vector<std::vector<T>>, 4>& x,
const std::array<std::vector<std::array<std::size_t, 2>>, 4>& shape)
{
std::array<std::vector<mdspan_t<const T, 2>>, 4> x1;
for (std::size_t i = 0; i < x.size(); ++i)
for (std::size_t j = 0; j < x[i].size(); ++j)
x1[i].push_back(mdspan_t<const T, 2>(x[i][j].data(), shape[i][j]));
return x1;
}
template <typename T>
std::array<std::vector<mdspan_t<const T, 4>>, 4>
to_mdspan(const std::array<std::vector<std::vector<T>>, 4>& M,
const std::array<std::vector<std::array<std::size_t, 4>>, 4>& shape)
{
std::array<std::vector<mdspan_t<const T, 4>>, 4> M1;
for (std::size_t i = 0; i < M.size(); ++i)
for (std::size_t j = 0; j < M[i].size(); ++j)
M1[i].push_back(mdspan_t<const T, 4>(M[i][j].data(), shape[i][j]));
return M1;
}
} // namespace impl
namespace element
{
template <typename T, std::size_t d>
using mdspan_t = impl::mdspan_t<T, d>;
template <std::floating_point T>
std::tuple<std::array<std::vector<std::vector<T>>, 4>,
std::array<std::vector<std::array<std::size_t, 2>>, 4>,
std::array<std::vector<std::vector<T>>, 4>,
std::array<std::vector<std::array<std::size_t, 4>>, 4>>
make_discontinuous(const std::array<std::vector<mdspan_t<const T, 2>>, 4>& x,
const std::array<std::vector<mdspan_t<const T, 4>>, 4>& M,
std::size_t tdim, std::size_t value_size);
} // namespace element
template <std::floating_point F>
class FiniteElement
{
template <typename T, std::size_t d>
using mdspan_t = md::mdspan<T, md::dextents<std::size_t, d>>;
public:
using scalar_type = F;
FiniteElement(element::family family, cell::type cell_type,
polyset::type poly_type, int degree,
const std::vector<std::size_t>& value_shape,
mdspan_t<const F, 2> wcoeffs,
const std::array<std::vector<mdspan_t<const F, 2>>, 4>& x,
const std::array<std::vector<mdspan_t<const F, 4>>, 4>& M,
int interpolation_nderivs, maps::type map_type,
sobolev::space sobolev_space, bool discontinuous,
int embedded_subdegree, int embedded_superdegree,
element::lagrange_variant lvariant,
element::dpc_variant dvariant,
std::vector<int> dof_ordering = {});
FiniteElement(const FiniteElement& element) = default;
FiniteElement(FiniteElement&& element) = default;
~FiniteElement() = default;
FiniteElement& operator=(const FiniteElement& element) = default;
FiniteElement& operator=(FiniteElement&& element) = default;
bool operator==(const FiniteElement& e) const;
std::size_t hash() const;
std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
std::size_t num_points) const
{
std::size_t ndsize = 1;
for (std::size_t i = 1; i <= nd; ++i)
ndsize *= (_cell_tdim + i);
for (std::size_t i = 1; i <= nd; ++i)
ndsize /= i;
std::size_t vs = std::accumulate(_value_shape.begin(), _value_shape.end(),
1, std::multiplies{});
std::size_t ndofs = _coeffs.second[0];
return {ndsize, num_points, ndofs, vs};
}
std::pair<std::vector<F>, std::array<std::size_t, 4>>
tabulate(int nd, impl::mdspan_t<const F, 2> x) const;
std::pair<std::vector<F>, std::array<std::size_t, 4>>
tabulate(int nd, std::span<const F> x,
std::array<std::size_t, 2> shape) const;
void tabulate(int nd, impl::mdspan_t<const F, 2> x,
mdspan_t<F, 4> basis) const;
void tabulate(int nd, std::span<const F> x, std::array<std::size_t, 2> xshape,
std::span<F> basis) const;
cell::type cell_type() const { return _cell_type; }
polyset::type polyset_type() const { return _poly_type; }
int degree() const { return _degree; }
int embedded_superdegree() const { return _embedded_superdegree; }
int embedded_subdegree() const { return _embedded_subdegree; }
const std::vector<std::size_t>& value_shape() const { return _value_shape; }
int dim() const { return _coeffs.second[0]; }
element::family family() const { return _family; }
element::lagrange_variant lagrange_variant() const
{
return _lagrange_variant;
}
element::dpc_variant dpc_variant() const { return _dpc_variant; }
maps::type map_type() const { return _map_type; }
sobolev::space sobolev_space() const { return _sobolev_space; }
bool discontinuous() const { return _discontinuous; }
bool dof_transformations_are_permutations() const
{
return _dof_transformations_are_permutations;
}
bool dof_transformations_are_identity() const
{
return _dof_transformations_are_identity;
}
std::pair<std::vector<F>, std::array<std::size_t, 3>>
push_forward(impl::mdspan_t<const F, 3> U, impl::mdspan_t<const F, 3> J,
std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
std::pair<std::vector<F>, std::array<std::size_t, 3>>
pull_back(impl::mdspan_t<const F, 3> u, impl::mdspan_t<const F, 3> J,
std::span<const F> detJ, impl::mdspan_t<const F, 3> K) const;
template <typename O, typename P, typename Q, typename R>
std::function<void(O&, const P&, const Q&, F, const R&)> map_fn() const
{
switch (_map_type)
{
case maps::type::identity:
return [](O& u, const P& U, const Q&, F, const R&)
{
assert(U.extent(0) == u.extent(0));
assert(U.extent(1) == u.extent(1));
for (std::size_t i = 0; i < U.extent(0); ++i)
for (std::size_t j = 0; j < U.extent(1); ++j)
u(i, j) = U(i, j);
};
case maps::type::covariantPiola:
return [](O& u, const P& U, const Q& J, F detJ, const R& K)
{ maps::covariant_piola(u, U, J, detJ, K); };
case maps::type::contravariantPiola:
return [](O& u, const P& U, const Q& J, F detJ, const R& K)
{ maps::contravariant_piola(u, U, J, detJ, K); };
case maps::type::doubleCovariantPiola:
return [](O& u, const P& U, const Q& J, F detJ, const R& K)
{ maps::double_covariant_piola(u, U, J, detJ, K); };
case maps::type::doubleContravariantPiola:
return [](O& u, const P& U, const Q& J, F detJ, const R& K)
{ maps::double_contravariant_piola(u, U, J, detJ, K); };
default:
throw std::runtime_error("Map not implemented");
}
}
const std::vector<std::vector<std::vector<int>>>& entity_dofs() const
{
return _edofs;
}
const std::vector<std::vector<std::vector<int>>>& entity_closure_dofs() const
{
return _e_closure_dofs;
}
std::pair<std::vector<F>, std::array<std::size_t, 3>>
base_transformations() const;
std::map<cell::type, std::pair<std::vector<F>, std::array<std::size_t, 3>>>
entity_transformations() const
{
return _entity_transformations;
}
void permute(std::span<std::int32_t> d, std::uint32_t cell_info) const
{
if (!_dof_transformations_are_permutations)
{
throw std::runtime_error(
"The DOF transformations for this element are not permutations");
}
if (_dof_transformations_are_identity)
return;
else
permute_data<std::int32_t, false>(d, 1, cell_info, _eperm);
}
void permute_inv(std::span<std::int32_t> d, std::uint32_t cell_info) const
{
if (!_dof_transformations_are_permutations)
{
throw std::runtime_error(
"The DOF transformations for this element are not permutations");
}
if (_dof_transformations_are_identity)
return;
else
permute_data<std::int32_t, true>(d, 1, cell_info, _eperm_inv);
}
void permute_subentity_closure(std::span<std::int32_t> d,
std::uint32_t cell_info,
cell::type entity_type, int entity_index) const
{
const int entity_dim = cell::topological_dimension(entity_type);
int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
std::uint32_t entity_info = 0;
switch (entity_dim)
{
case 0:
return;
case 1:
entity_info = cell_info >> (face_start + entity_index) & 1;
break;
case 2:
entity_info = cell_info >> (3 * entity_index) & 7;
break;
default:
throw std::runtime_error("Unsupported cell dimension");
}
permute_subentity_closure(d, entity_info, entity_type);
}
void permute_subentity_closure_inv(std::span<std::int32_t> d,
std::uint32_t cell_info,
cell::type entity_type,
int entity_index) const
{
const int entity_dim = cell::topological_dimension(entity_type);
int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
std::uint32_t entity_info;
switch (entity_dim)
{
case 0:
return;
case 1:
entity_info = cell_info >> (face_start + entity_index) & 1;
break;
case 2:
entity_info = cell_info >> (3 * entity_index) & 7;
break;
default:
throw std::runtime_error("Unsupported cell dimension");
}
permute_subentity_closure_inv(d, entity_info, entity_type);
}
void permute_subentity_closure(std::span<std::int32_t> d,
std::uint32_t entity_info,
cell::type entity_type) const
{
if (!_dof_transformations_are_permutations)
{
throw std::runtime_error(
"The DOF transformations for this element are not permutations");
}
const int entity_dim = cell::topological_dimension(entity_type);
if (entity_dim == 0)
return;
auto& perm = _subentity_closure_perm.at(entity_type);
if (entity_dim == 1)
{
if (entity_info & 1)
{
precompute::apply_permutation(perm[0], d);
}
}
else if (entity_dim == 2)
{
// Rotate a face
for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
{
precompute::apply_permutation(perm[0], d);
}
// Reflect a face (post rotate)
if (entity_info & 1)
{
precompute::apply_permutation(perm[1], d);
}
}
else
{
throw std::runtime_error(
"Invalid dimension for permute_subentity_closure");
}
}
void permute_subentity_closure_inv(std::span<std::int32_t> d,
std::uint32_t entity_info,
cell::type entity_type) const
{
if (!_dof_transformations_are_permutations)
{
throw std::runtime_error(
"The DOF transformations for this element are not permutations");
}
const int entity_dim = cell::topological_dimension(entity_type);
if (entity_dim == 0)
return;
auto& perm = _subentity_closure_perm_inv.at(entity_type);
if (entity_dim == 1)
{
if (entity_info & 1)
{
precompute::apply_permutation(perm[0], d);
}
}
else if (entity_dim == 2)
{
// Reflect a face (pre rotate)
if (entity_info & 1)
{
precompute::apply_permutation(perm[1], d);
}
// Rotate a face
for (std::uint32_t r = 0; r < (entity_info >> 1 & 3); ++r)
{
precompute::apply_permutation(perm[0], d);
}
}
else
{
throw std::runtime_error(
"Invalid dimension for permute_subentity_closure");
}
}
template <typename T>
void T_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tt_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tt_inv_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tinv_apply(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tt_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void T_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tinv_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
template <typename T>
void Tt_inv_apply_right(std::span<T> u, int n, std::uint32_t cell_info) const;
const std::pair<std::vector<F>, std::array<std::size_t, 2>>& points() const
{
return _points;
}
const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
interpolation_matrix() const
{
return _matM;
}
const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
dual_matrix() const
{
return _dual_matrix;
}
const std::pair<std::vector<F>, std::array<std::size_t, 2>>& wcoeffs() const
{
return _wcoeffs;
}
const std::array<
std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>, 4>&
x() const
{
return _x;
}
const std::array<
std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>, 4>&
M() const
{
return _M;
}
const std::pair<std::vector<F>, std::array<std::size_t, 2>>&
coefficient_matrix() const
{
return _coeffs;
}
bool has_tensor_product_factorisation() const
{
return !_tensor_factors.empty();
}
std::vector<std::vector<FiniteElement<F>>>
get_tensor_product_representation() const
{
if (!has_tensor_product_factorisation())
throw std::runtime_error("Element has no tensor product representation.");
return _tensor_factors;
}
bool interpolation_is_identity() const { return _interpolation_is_identity; }
int interpolation_nderivs() const { return _interpolation_nderivs; }
const std::vector<int>& dof_ordering() const { return _dof_ordering; }
private:
template <typename T, bool post>
void permute_data(
std::span<T> data, int block_size, std::uint32_t cell_info,
const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
const;
using array2_t = std::pair<std::vector<F>, std::array<std::size_t, 2>>;
using array3_t = std::pair<std::vector<F>, std::array<std::size_t, 3>>;
using trans_data_t
= std::vector<std::pair<std::vector<std::size_t>, array2_t>>;
template <typename T, bool post, typename OP>
void
transform_data(std::span<T> data, int block_size, std::uint32_t cell_info,
const std::map<cell::type, trans_data_t>& etrans, OP op) const;
// Cell type
cell::type _cell_type;
// Polyset type
polyset::type _poly_type;
// Topological dimension of the cell
std::size_t _cell_tdim;
// Topological dimension of the cell
std::vector<std::vector<cell::type>> _cell_subentity_types;
// Finite element family
element::family _family;
// Lagrange variant
element::lagrange_variant _lagrange_variant;
// DPC variant
element::dpc_variant _dpc_variant;
// Degree that was input when creating the element
int _degree;
// Degree
int _interpolation_nderivs;
// Highest degree polynomial in element's polyset
int _embedded_superdegree;
// Highest degree space that is a subspace of element's polyset
int _embedded_subdegree;
// Value shape
std::vector<std::size_t> _value_shape;
maps::type _map_type;
sobolev::space _sobolev_space;
// Shape function coefficient of expansion sets on cell. If shape
// function is given by @f$\psi_i = \sum_{k} \phi_{k}
// \alpha^{i}_{k}@f$, then _coeffs(i, j) = @f$\alpha^i_k@f$. ie
// _coeffs.row(i) are the expansion coefficients for shape function i
// (@f$\psi_{i}@f$).
std::pair<std::vector<F>, std::array<std::size_t, 2>> _coeffs;
// Dofs associated with each cell (sub-)entity
std::vector<std::vector<std::vector<int>>> _edofs;
// Dofs associated with the closdure of each cell (sub-)entity
std::vector<std::vector<std::vector<int>>> _e_closure_dofs;
// Entity transformations
std::map<cell::type, array3_t> _entity_transformations;
// Set of points used for point evaluation
// Experimental - currently used for an implementation of
// "tabulate_dof_coordinates" Most useful for Lagrange. This may change or go
// away. For non-Lagrange elements, these points will be used in combination
// with _interpolation_matrix to perform interpolation
std::pair<std::vector<F>, std::array<std::size_t, 2>> _points;
// Interpolation points on the cell. The shape is (entity_dim, num
// entities of given dimension, num_points, tdim)
std::array<std::vector<std::pair<std::vector<F>, std::array<std::size_t, 2>>>,
4>
_x;
std::pair<std::vector<F>, std::array<std::size_t, 2>> _matM;
// Indicates whether or not the DOF transformations are all
// permutations
bool _dof_transformations_are_permutations;
// Indicates whether or not the DOF transformations are all identity
bool _dof_transformations_are_identity;
// The entity permutations (factorised). This will only be set if
// _dof_transformations_are_permutations is True and
// _dof_transformations_are_identity is False
std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm;
// The reverse entity permutations (factorised). This will only be set
// if _dof_transformations_are_permutations is True and
// _dof_transformations_are_identity is False
std::map<cell::type, std::vector<std::vector<std::size_t>>> _eperm_inv;
// The entity transformations in precomputed form
std::map<cell::type, trans_data_t> _etrans;
// The transposed entity transformations in precomputed form
std::map<cell::type, trans_data_t> _etransT;
// The inverse entity transformations in precomputed form
std::map<cell::type, trans_data_t> _etrans_inv;
// The inverse transpose entity transformations in precomputed form
std::map<cell::type, trans_data_t> _etrans_invT;
// The subentity closure permutations (factorised). This will only be set if
// _dof_transformations_are_permutations is True
std::map<cell::type, std::vector<std::vector<std::size_t>>>
_subentity_closure_perm;
// The inverse subentity closure permutations (factorised). This will only be
// set if _dof_transformations_are_permutations is True
std::map<cell::type, std::vector<std::vector<std::size_t>>>
_subentity_closure_perm_inv;
// Indicates whether or not this is the discontinuous version of the
// element
bool _discontinuous;
// The dual matrix
std::pair<std::vector<F>, std::array<std::size_t, 2>> _dual_matrix;
// Dof reordering for different element dof layout compatibility.
// The reference basix layout is ordered by entity, i.e. dofs on
// vertices, followed by edges, faces, then internal dofs.
// _dof_ordering stores the map to the new order required, e.g.
// for a P2 triangle, _dof_ordering=[0 3 5 1 2 4] will place
// dofs 0, 3, 5 on the vertices and 1, 2, 4, on the edges.
std::vector<int> _dof_ordering;
// Tensor product representation
// Entries of tuple are (list of elements on an interval, permutation
// of DOF numbers)
// @todo: For vector-valued elements, a tensor product type and a
// scaling factor may additionally be needed.
std::vector<std::vector<FiniteElement>> _tensor_factors;
// Is the interpolation matrix an identity?
bool _interpolation_is_identity;
// The coefficients that define the polynomial set in terms of the
// orthonormal polynomials
std::pair<std::vector<F>, std::array<std::size_t, 2>> _wcoeffs;
// Interpolation matrices for each entity
using array4_t
= std::vector<std::pair<std::vector<F>, std::array<std::size_t, 4>>>;
std::array<array4_t, 4> _M;
};
template <std::floating_point T>
FiniteElement<T> create_custom_element(
cell::type cell_type, const std::vector<std::size_t>& value_shape,
impl::mdspan_t<const T, 2> wcoeffs,
const std::array<std::vector<impl::mdspan_t<const T, 2>>, 4>& x,
const std::array<std::vector<impl::mdspan_t<const T, 4>>, 4>& M,
int interpolation_nderivs, maps::type map_type,
sobolev::space sobolev_space, bool discontinuous, int embedded_subdegree,
int embedded_superdegree, polyset::type poly_type);
template <std::floating_point T>
FiniteElement<T> create_element(element::family family, cell::type cell,
int degree, element::lagrange_variant lvariant,
element::dpc_variant dvariant,
bool discontinuous,
std::vector<int> dof_ordering = {});
std::vector<int> tp_dof_ordering(element::family family, cell::type cell,
int degree, element::lagrange_variant lvariant,
element::dpc_variant dvariant,
bool discontinuous);
std::vector<int> lex_dof_ordering(element::family family, cell::type cell,
int degree, element::lagrange_variant lvariant,
element::dpc_variant dvariant,
bool discontinuous);
template <std::floating_point T>
std::vector<std::vector<FiniteElement<T>>>
tp_factors(element::family family, cell::type cell, int degree,
element::lagrange_variant lvariant, element::dpc_variant dvariant,
bool discontinuous, const std::vector<int>& dof_ordering);
template <std::floating_point T>
FiniteElement<T>
create_tp_element(element::family family, cell::type cell, int degree,
element::lagrange_variant lvariant,
element::dpc_variant dvariant, bool discontinuous);
std::string version();
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T, bool post>
void FiniteElement<F>::permute_data(
std::span<T> data, int block_size, std::uint32_t cell_info,
const std::map<cell::type, std::vector<std::vector<std::size_t>>>& eperm)
const
{
if (_cell_tdim >= 2)
{
// This assumes 3 bits are used per face. This will need updating if 3D
// cells with faces with more than 4 sides are implemented
int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
// Permute DOFs on edges
{
auto& trans = eperm.at(cell::type::interval)[0];
for (std::size_t e = 0; e < _edofs[1].size(); ++e)
{
// Reverse an edge
if (cell_info >> (face_start + e) & 1)
{
precompute::apply_permutation_mapped(trans, data, _edofs[1][e],
block_size);
}
}
}
if (_cell_tdim == 3)
{
// Permute DOFs on faces
for (std::size_t f = 0; f < _edofs[2].size(); ++f)
{
auto& trans = eperm.at(_cell_subentity_types[2][f]);
// Reflect a face (pre rotate)
if (!post and cell_info >> (3 * f) & 1)
{
precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
block_size);
}
// Rotate a face
for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
{
precompute::apply_permutation_mapped(trans[0], data, _edofs[2][f],
block_size);
}
// Reflect a face (post rotate)
if (post and cell_info >> (3 * f) & 1)
{
precompute::apply_permutation_mapped(trans[1], data, _edofs[2][f],
block_size);
}
}
}
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T, bool post, typename OP>
void FiniteElement<F>::transform_data(
std::span<T> data, int block_size, std::uint32_t cell_info,
const std::map<cell::type, trans_data_t>& etrans, OP op) const
{
if (_cell_tdim >= 2)
{
// This assumes 3 bits are used per face. This will need updating if
// 3D cells with faces with more than 4 sides are implemented
int face_start = _cell_tdim == 3 ? 3 * _edofs[2].size() : 0;
int dofstart = 0;
for (auto& edofs0 : _edofs[0])
dofstart += edofs0.size();
// Transform DOFs on edges
{
auto& [v_size_t, matrix] = etrans.at(cell::type::interval)[0];
for (std::size_t e = 0; e < _edofs[1].size(); ++e)
{
// Reverse an edge
if (cell_info >> (face_start + e) & 1)
{
op(std::span(v_size_t),
mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
dofstart, block_size);
}
dofstart += _edofs[1][e].size();
}
}
if (_cell_tdim == 3)
{
// Permute DOFs on faces
for (std::size_t f = 0; f < _edofs[2].size(); ++f)
{
auto& trans = etrans.at(_cell_subentity_types[2][f]);
// Reflect a face (pre rotation)
if (!post and cell_info >> (3 * f) & 1)
{
const auto& m = trans[1];
const auto& v_size_t = std::get<0>(m);
const auto& matrix = std::get<1>(m);
op(std::span(v_size_t),
mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
dofstart, block_size);
}
// Rotate a face
for (std::uint32_t r = 0; r < (cell_info >> (3 * f + 1) & 3); ++r)
{
const auto& m = trans[0];
const auto& v_size_t = std::get<0>(m);
const auto& matrix = std::get<1>(m);
op(std::span(v_size_t),
mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
dofstart, block_size);
}
// Reflect a face (post rotation)
if (post and cell_info >> (3 * f) & 1)
{
const auto& m = trans[1];
const auto& v_size_t = std::get<0>(m);
const auto& matrix = std::get<1>(m);
op(std::span(v_size_t),
mdspan_t<const F, 2>(matrix.first.data(), matrix.second), data,
dofstart, block_size);
}
dofstart += _edofs[2][f].size();
}
}
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::T_apply(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
if (_dof_transformations_are_permutations)
permute_data<T, false>(u, n, cell_info, _eperm);
else
{
transform_data<T, false>(u, n, cell_info, _etrans,
precompute::apply_matrix<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tt_apply(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
permute_data<T, true>(u, n, cell_info, _eperm_inv);
else
{
transform_data<T, true>(u, n, cell_info, _etransT,
precompute::apply_matrix<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tt_inv_apply(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
permute_data<T, false>(u, n, cell_info, _eperm);
else
{
transform_data<T, false>(u, n, cell_info, _etrans_invT,
precompute::apply_matrix<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tinv_apply(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
permute_data<T, true>(u, n, cell_info, _eperm_inv);
else
{
transform_data<T, true>(u, n, cell_info, _etrans_inv,
precompute::apply_matrix<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tt_apply_right(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
{
assert(u.size() % n == 0);
const int step = u.size() / n;
for (int i = 0; i < n; ++i)
{
std::span<T> dblock(u.data() + i * step, step);
permute_data<T, false>(dblock, 1, cell_info, _eperm);
}
}
else
{
transform_data<T, false>(u, n, cell_info, _etrans,
precompute::apply_tranpose_matrix_right<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tinv_apply_right(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
{
assert(u.size() % n == 0);
const int step = u.size() / n;
for (int i = 0; i < n; ++i)
{
std::span<T> dblock(u.data() + i * step, step);
permute_data<T, false>(dblock, 1, cell_info, _eperm);
}
}
else
{
transform_data<T, false>(u, n, cell_info, _etrans_invT,
precompute::apply_tranpose_matrix_right<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::T_apply_right(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
{
assert(u.size() % n == 0);
const int step = u.size() / n;
for (int i = 0; i < n; ++i)
{
std::span<T> dblock(u.data() + i * step, step);
permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
}
}
else
{
transform_data<T, true>(u, n, cell_info, _etransT,
precompute::apply_tranpose_matrix_right<F, T>);
}
}
//-----------------------------------------------------------------------------
template <std::floating_point F>
template <typename T>
void FiniteElement<F>::Tt_inv_apply_right(std::span<T> u, int n,
std::uint32_t cell_info) const
{
if (_dof_transformations_are_identity)
return;
else if (_dof_transformations_are_permutations)
{
assert(u.size() % n == 0);
const int step = u.size() / n;
for (int i = 0; i < n; ++i)
{
std::span<T> dblock(u.data() + i * step, step);
permute_data<T, true>(dblock, 1, cell_info, _eperm_inv);
}
}
else
{
transform_data<T, true>(u, n, cell_info, _etrans_inv,
precompute::apply_tranpose_matrix_right<F, T>);
}
}
//-----------------------------------------------------------------------------
} // namespace basix
basix::FiniteElement
A finite element.
Definition finite-element.h:137
basix::FiniteElement::FiniteElement
FiniteElement(element::family family, cell::type cell_type, polyset::type poly_type, int degree, const std::vector< std::size_t > &value_shape, mdspan_t< const F, 2 > wcoeffs, const std::array< std::vector< mdspan_t< const F, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const F, 4 > >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int embedded_subdegree, int embedded_superdegree, element::lagrange_variant lvariant, element::dpc_variant dvariant, std::vector< int > dof_ordering={})
Construct a finite element.
basix::FiniteElement::M
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 4 > > >, 4 > & M() const
Get the interpolation matrices for each subentity.
Definition finite-element.h:1325
basix::FiniteElement::Tinv_apply_right
void Tinv_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the inverse of the operator applied by T_apply().
Definition finite-element.h:1914
basix::FiniteElement::base_transformations
std::pair< std::vector< F >, std::array< std::size_t, 3 > > base_transformations() const
Get the base transformations.
Definition finite-element.cpp:1888
basix::FiniteElement::T_apply
void T_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Transform basis functions from the reference element ordering and orientation to the globally consist...
Definition finite-element.h:1825
basix::FiniteElement::Tt_apply
void Tt_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the transpose of the operator applied by T_apply().
Definition finite-element.h:1842
basix::FiniteElement::permute_subentity_closure_inv
void permute_subentity_closure_inv(std::span< std::int32_t > d, std::uint32_t cell_info, cell::type entity_type, int entity_index) const
Perform the inverse of the operation applied by permute_subentity_closure().
Definition finite-element.h:892
basix::FiniteElement::dof_transformations_are_identity
bool dof_transformations_are_identity() const
Indicates is the dof transformations are all the identity.
Definition finite-element.h:557
basix::FiniteElement::operator==
bool operator==(const FiniteElement &e) const
Check if two elements are the same.
Definition finite-element.cpp:1710
basix::FiniteElement::entity_transformations
std::map< cell::type, std::pair< std::vector< F >, std::array< std::size_t, 3 > > > entity_transformations() const
Return the entity dof transformation matrices.
Definition finite-element.h:768
basix::FiniteElement::dof_ordering
const std::vector< int > & dof_ordering() const
Get dof layout.
Definition finite-element.h:1384
basix::FiniteElement::permute_subentity_closure_inv
void permute_subentity_closure_inv(std::span< std::int32_t > d, std::uint32_t entity_info, cell::type entity_type) const
Perform the inverse of the operation applied by permute_subentity_closure().
Definition finite-element.h:988
basix::FiniteElement::embedded_subdegree
int embedded_subdegree() const
Highest degree n such that a Lagrange (or vector Lagrange) element of degree n is a subspace of this ...
Definition finite-element.h:505
basix::FiniteElement::entity_dofs
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const
Get the dofs on each topological entity: (vertices, edges, faces, cell) in that order.
Definition finite-element.h:660
basix::FiniteElement::T_apply_right
void T_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the operator applied by T_apply().
Definition finite-element.h:1938
basix::FiniteElement::FiniteElement
FiniteElement(FiniteElement &&element)=default
Move constructor.
basix::FiniteElement::has_tensor_product_factorisation
bool has_tensor_product_factorisation() const
Indicates whether or not this element can be represented as a product of elements defined on lower-di...
Definition finite-element.h:1349
basix::FiniteElement::interpolation_nderivs
int interpolation_nderivs() const
The number of derivatives needed when interpolating.
Definition finite-element.h:1381
basix::FiniteElement::tabulate
std::pair< std::vector< F >, std::array< std::size_t, 4 > > tabulate(int nd, impl::mdspan_t< const F, 2 > x) const
Compute basis values and derivatives at set of points.
Definition finite-element.cpp:1796
basix::FiniteElement::embedded_superdegree
int embedded_superdegree() const
Lowest degree n such that the highest degree polynomial in this element is contained in a Lagrange (o...
Definition finite-element.h:500
basix::FiniteElement::dim
int dim() const
Dimension of the finite element space.
Definition finite-element.h:518
basix::FiniteElement::permute
void permute(std::span< std::int32_t > d, std::uint32_t cell_info) const
Permute indices associated with degree-of-freedoms on the reference element ordering to the globally ...
Definition finite-element.h:795
basix::FiniteElement::hash
std::size_t hash() const
Get a unique hash of this element.
Definition finite-element.cpp:1749
basix::FiniteElement::scalar_type
F scalar_type
Scalar type.
Definition finite-element.h:143
basix::FiniteElement::permute_inv
void permute_inv(std::span< std::int32_t > d, std::uint32_t cell_info) const
Perform the inverse of the operation applied by permute().
Definition finite-element.h:826
basix::FiniteElement::Tt_apply_right
void Tt_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the transpose of the operator applied by T_apply().
Definition finite-element.h:1890
basix::FiniteElement::value_shape
const std::vector< std::size_t > & value_shape() const
Element value tensor shape.
Definition finite-element.h:512
basix::FiniteElement::Tt_inv_apply
void Tt_inv_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the inverse transpose of the operator applied by T_apply().
Definition finite-element.h:1858
basix::FiniteElement::permute_subentity_closure
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t cell_info, cell::type entity_type, int entity_index) const
Permute indices associated with degree-of-freedoms on the closure of a sub-entity of the reference el...
Definition finite-element.h:855
basix::FiniteElement::lagrange_variant
element::lagrange_variant lagrange_variant() const
Lagrange variant of the element.
Definition finite-element.h:526
basix::FiniteElement::Tt_inv_apply_right
void Tt_inv_apply_right(std::span< T > u, int n, std::uint32_t cell_info) const
Right(post)-apply the transpose inverse of the operator applied by T_apply().
Definition finite-element.h:1962
basix::FiniteElement::x
const std::array< std::vector< std::pair< std::vector< F >, std::array< std::size_t, 2 > > >, 4 > & x() const
Get the interpolation points for each subentity.
Definition finite-element.h:1282
basix::FiniteElement::degree
int degree() const
Get the element polynomial degree.
Definition finite-element.h:494
basix::FiniteElement::permute_subentity_closure
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t entity_info, cell::type entity_type) const
Permute indices associated with degree-of-freedoms on the closure of a sub-entity of the reference el...
Definition finite-element.h:933
basix::FiniteElement::Tinv_apply
void Tinv_apply(std::span< T > u, int n, std::uint32_t cell_info) const
Apply the inverse of the operator applied by T_apply().
Definition finite-element.h:1874
basix::FiniteElement::points
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & points() const
Return the interpolation points.
Definition finite-element.h:1162
basix::FiniteElement::sobolev_space
sobolev::space sobolev_space() const
Underlying Sobolev space for this element.
Definition finite-element.h:541
basix::FiniteElement::entity_closure_dofs
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const
Get the dofs on the closure of each topological entity: (vertices, edges, faces, cell) in that order.
Definition finite-element.h:674
basix::FiniteElement::FiniteElement
FiniteElement(const FiniteElement &element)=default
Copy constructor.
basix::FiniteElement::tabulate_shape
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Array shape for tabulate basis values and derivatives at set of points.
Definition finite-element.h:363
basix::FiniteElement::operator=
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment operator.
basix::FiniteElement::interpolation_matrix
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & interpolation_matrix() const
Return a matrix of weights interpolation.
Definition finite-element.h:1219
basix::FiniteElement::polyset_type
polyset::type polyset_type() const
Get the element polyset type.
Definition finite-element.h:490
basix::FiniteElement::wcoeffs
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & wcoeffs() const
Get the coefficients that define the polynomial set in terms of the orthonormal polynomials.
Definition finite-element.h:1271
basix::FiniteElement::pull_back
std::pair< std::vector< F >, std::array< std::size_t, 3 > > pull_back(impl::mdspan_t< const F, 3 > u, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Map function values from a physical cell to the reference.
Definition finite-element.cpp:1993
basix::FiniteElement::dof_transformations_are_permutations
bool dof_transformations_are_permutations() const
Indicates if the degree-of-freedom transformations are all permutations.
Definition finite-element.h:551
basix::FiniteElement::push_forward
std::pair< std::vector< F >, std::array< std::size_t, 3 > > push_forward(impl::mdspan_t< const F, 3 > U, impl::mdspan_t< const F, 3 > J, std::span< const F > detJ, impl::mdspan_t< const F, 3 > K) const
Map function values from the reference to a physical cell.
Definition finite-element.cpp:1957
basix::FiniteElement::cell_type
cell::type cell_type() const
Get the element cell type.
Definition finite-element.h:486
basix::FiniteElement::interpolation_is_identity
bool interpolation_is_identity() const
Indicates whether or not the interpolation matrix for this element is an identity matrix.
Definition finite-element.h:1378
basix::FiniteElement::discontinuous
bool discontinuous() const
Indicates whether this element is the discontinuous variant.
Definition finite-element.h:547
basix::FiniteElement::get_tensor_product_representation
std::vector< std::vector< FiniteElement< F > > > get_tensor_product_representation() const
Get the tensor product representation of this element.
Definition finite-element.h:1367
basix::FiniteElement::map_type
maps::type map_type() const
Map type for the element.
Definition finite-element.h:537
basix::FiniteElement::dual_matrix
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & dual_matrix() const
Get the dual matrix.
Definition finite-element.h:1230
basix::FiniteElement::dpc_variant
element::dpc_variant dpc_variant() const
DPC variant of the element.
Definition finite-element.h:533
basix::FiniteElement::coefficient_matrix
const std::pair< std::vector< F >, std::array< std::size_t, 2 > > & coefficient_matrix() const
Get the matrix of coefficients.
Definition finite-element.h:1334
basix::FiniteElement::operator=
FiniteElement & operator=(const FiniteElement &element)=default
Assignment operator.
basix::FiniteElement::family
element::family family() const
The finite element family.
Definition finite-element.h:522
basix::FiniteElement::~FiniteElement
~FiniteElement()=default
Destructor.
basix::FiniteElement::map_fn
std::function< void(O &, const P &, const Q &, F, const R &)> map_fn() const
Return a function that performs the appropriate push-forward/pull-back for the element type.
Definition finite-element.h:623
basix::cell
Information about reference cells.
Definition cell.h:17
basix::cell::type
type
Cell type.
Definition cell.h:21
basix::cell::topological_dimension
int topological_dimension(cell::type celltype)
Definition cell.cpp:294
basix::element
Interfaces for creating finite elements.
Definition e-brezzi-douglas-marini.h:13
basix::element::mdspan_t
impl::mdspan_t< T, d > mdspan_t
Typedef for mdspan.
Definition finite-element.h:102
basix::element::make_discontinuous
std::tuple< std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 2 > >, 4 >, std::array< std::vector< std::vector< T > >, 4 >, std::array< std::vector< std::array< std::size_t, 4 > >, 4 > > make_discontinuous(const std::array< std::vector< mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< mdspan_t< const T, 4 > >, 4 > &M, std::size_t tdim, std::size_t value_size)
Definition finite-element.cpp:819
basix::element::lagrange_variant
lagrange_variant
Variants of a Lagrange space that can be created.
Definition element-families.h:12
basix::element::dpc_variant
dpc_variant
Definition element-families.h:32
basix::element::family
family
Available element families.
Definition element-families.h:45
basix::maps::covariant_piola
void covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Covariant Piola map.
Definition maps.h:62
basix::maps::contravariant_piola
void contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Contravariant Piola map.
Definition maps.h:82
basix::maps::double_contravariant_piola
void double_contravariant_piola(O &&r, const P &U, const Q &J, double detJ, const R &)
Double contravariant Piola map.
Definition maps.h:132
basix::maps::double_covariant_piola
void double_covariant_piola(O &&r, const P &U, const Q &, double, const R &K)
Double covariant Piola map.
Definition maps.h:104
basix::maps::type
type
Map type.
Definition maps.h:40
basix::polyset::type
type
Cell type.
Definition polyset.h:137
basix::precompute::apply_permutation
void apply_permutation(std::span< const std::size_t > perm, std::span< E > data, std::size_t offset=0, std::size_t n=1)
Apply a (precomputed) permutation .
Definition precompute.h:137
basix::precompute::apply_permutation_mapped
void apply_permutation_mapped(std::span< const std::size_t > perm, std::span< E > data, std::span< const int > emap, std::size_t n=1)
Permutation of mapped data.
Definition precompute.h:147
basix::precompute::apply_tranpose_matrix_right
void apply_tranpose_matrix_right(std::span< const std::size_t > v_size_t, md::mdspan< const T, md::dextents< std::size_t, 2 > > M, std::span< E > data, std::size_t offset=0, std::size_t n=1)
Apply a (precomputed) matrix to some transposed data.
Definition precompute.h:276
basix::precompute::apply_matrix
void apply_matrix(std::span< const std::size_t > v_size_t, md::mdspan< const T, md::dextents< std::size_t, 2 > > M, std::span< E > data, std::size_t offset=0, std::size_t n=1)
Apply a (precomputed) matrix.
Definition precompute.h:236
basix::sobolev::space
space
Sobolev space type.
Definition sobolev-spaces.h:13
basix
Basix: FEniCS runtime basis evaluation library.
Definition cell.h:17
basix::create_element
FiniteElement< T > create_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous, std::vector< int > dof_ordering={})
Definition finite-element.cpp:190
basix::tp_factors
std::vector< std::vector< FiniteElement< T > > > tp_factors(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous, const std::vector< int > &dof_ordering)
Definition finite-element.cpp:346
basix::tp_dof_ordering
std::vector< int > tp_dof_ordering(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous)
Definition finite-element.cpp:396
basix::version
std::string version()
Definition finite-element.cpp:2027
basix::lex_dof_ordering
std::vector< int > lex_dof_ordering(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous)
Definition finite-element.cpp:512
basix::create_custom_element
FiniteElement< T > create_custom_element(cell::type cell_type, const std::vector< std::size_t > &value_shape, impl::mdspan_t< const T, 2 > wcoeffs, const std::array< std::vector< impl::mdspan_t< const T, 2 > >, 4 > &x, const std::array< std::vector< impl::mdspan_t< const T, 4 > >, 4 > &M, int interpolation_nderivs, maps::type map_type, sobolev::space sobolev_space, bool discontinuous, int embedded_subdegree, int embedded_superdegree, polyset::type poly_type)
Definition finite-element.cpp:904
basix::create_tp_element
FiniteElement< T > create_tp_element(element::family family, cell::type cell, int degree, element::lagrange_variant lvariant, element::dpc_variant dvariant, bool discontinuous)
Definition finite-element.cpp:327