DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
FiniteElement.h
1// Copyright (C) 2020-2024 Garth N. Wells and Matthew W. Scroggs
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "traits.h"
10#include <array>
11#include <basix/finite-element.h>
12#include <concepts>
13#include <cstdint>
14#include <dolfinx/mesh/cell_types.h>
15#include <functional>
16#include <memory>
17#include <optional>
18#include <span>
19#include <utility>
20#include <vector>
21
22namespace dolfinx::fem
23{
32
35template <std::floating_point T>
37{
38 std::reference_wrapper<const basix::FiniteElement<T>>
40 std::optional<std::vector<std::size_t>> value_shape
41 = std::nullopt;
42 bool symmetry = false;
44};
45
47template <typename U, typename V, typename W>
48BasixElementData(U element, V bs, W symmetry)
50
55template <std::floating_point T>
57{
58public:
60 using geometry_type = T;
61
71 std::optional<std::vector<std::size_t>> value_shape
72 = std::nullopt,
73 bool symmetric = false);
74
84
103 const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
104 elements);
105
112 FiniteElement(mesh::CellType cell_type, std::span<const geometry_type> points,
113 std::array<std::size_t, 2> pshape,
114 std::vector<std::size_t> value_shape = {},
115 bool symmetric = false);
116
118 FiniteElement(const FiniteElement& element) = delete;
119
121 FiniteElement(FiniteElement&& element) = default;
122
124 ~FiniteElement() = default;
125
127 FiniteElement& operator=(const FiniteElement& element) = delete;
128
130 FiniteElement& operator=(FiniteElement&& element) = default;
131
136 bool operator==(const FiniteElement& e) const;
137
142 bool operator!=(const FiniteElement& e) const;
143
145 mesh::CellType cell_type() const noexcept;
146
153 std::string signature() const noexcept;
154
162 int space_dimension() const noexcept;
163
174 int block_size() const noexcept;
175
195 int value_size() const;
196
206 std::span<const std::size_t> value_shape() const;
207
223 int reference_value_size() const;
224
236 std::span<const std::size_t> reference_value_shape() const;
237
239 const std::vector<std::vector<std::vector<int>>>&
240 entity_dofs() const noexcept;
241
244 const std::vector<std::vector<std::vector<int>>>&
245 entity_closure_dofs() const noexcept;
246
248 bool symmetric() const;
249
262 void tabulate(std::span<geometry_type> values,
263 std::span<const geometry_type> X,
264 std::array<std::size_t, 2> shape, int order) const;
265
276 std::pair<std::vector<geometry_type>, std::array<std::size_t, 4>>
277 tabulate(std::span<const geometry_type> X, std::array<std::size_t, 2> shape,
278 int order) const;
279
282 int num_sub_elements() const noexcept;
283
291 bool is_mixed() const noexcept;
292
294 const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
295 sub_elements() const noexcept;
296
298 std::shared_ptr<const FiniteElement<geometry_type>>
299 extract_sub_element(const std::vector<int>& component) const;
300
303 const basix::FiniteElement<geometry_type>& basix_element() const;
304
306 basix::maps::type map_type() const;
307
314 bool interpolation_ident() const noexcept;
315
320 bool map_ident() const noexcept;
321
332 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
333 interpolation_points() const;
334
344 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
346
359 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
361
375 bool needs_dof_transformations() const noexcept;
376
391 bool needs_dof_permutations() const noexcept;
392
432 template <typename U>
433 std::function<void(std::span<U>, std::span<const std::uint32_t>, std::int32_t,
434 int)>
435 dof_transformation_fn(doftransform ttype, bool scalar_element = false) const
436 {
438 {
439 // If no permutation needed, return function that does nothing
440 return [](std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)
441 {
442 // Do nothing
443 };
444 }
445
446 if (!_sub_elements.empty())
447 {
448 if (!_reference_value_shape) // Mixed element
449 {
450 std::vector<std::function<void(
451 std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
452 sub_element_fns;
453 std::vector<int> dims;
454 for (std::size_t i = 0; i < _sub_elements.size(); ++i)
455 {
456 sub_element_fns.push_back(
457 _sub_elements[i]->template dof_transformation_fn<U>(ttype));
458 dims.push_back(_sub_elements[i]->space_dimension());
459 }
460
461 return [dims, sub_element_fns](std::span<U> data,
462 std::span<const std::uint32_t> cell_info,
463 std::int32_t cell, int block_size)
464 {
465 std::size_t offset = 0;
466 for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
467 {
468 const std::size_t width = dims[e] * block_size;
469 sub_element_fns[e](data.subspan(offset, width), cell_info, cell,
470 block_size);
471 offset += width;
472 }
473 };
474 }
475 else if (!scalar_element)
476 {
477 // Blocked element
478 std::function<void(std::span<U>, std::span<const std::uint32_t>,
479 std::int32_t, int)>
480 sub_fn
481 = _sub_elements.front()->template dof_transformation_fn<U>(ttype);
482 const int ebs = _bs;
483 return [ebs, sub_fn](std::span<U> data,
484 std::span<const std::uint32_t> cell_info,
485 std::int32_t cell, int data_block_size)
486 { sub_fn(data, cell_info, cell, ebs * data_block_size); };
487 }
488 }
489
490 switch (ttype)
491 {
493 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
494 std::int32_t cell, int block_size)
495 { Tt_inv_apply(data, cell_info[cell], block_size); };
497 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
498 std::int32_t cell, int block_size)
499 { Tt_apply(data, cell_info[cell], block_size); };
501 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
502 std::int32_t cell, int block_size)
503 { Tinv_apply(data, cell_info[cell], block_size); };
505 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
506 std::int32_t cell, int block_size)
507 { T_apply(data, cell_info[cell], block_size); };
508 default:
509 throw std::runtime_error("Unknown transformation type");
510 }
511 }
512
536 template <typename U>
537 std::function<void(std::span<U>, std::span<const std::uint32_t>, std::int32_t,
538 int)>
540 bool scalar_element = false) const
541 {
543 {
544 // If no permutation needed, return function that does nothing
545 return [](std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)
546 {
547 // Do nothing
548 };
549 }
550 else if (!_sub_elements.empty())
551 {
552 if (!_reference_value_shape) // Mixed element
553 {
554 std::vector<std::function<void(
555 std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
556 sub_element_fns;
557 for (std::size_t i = 0; i < _sub_elements.size(); ++i)
558 {
559 sub_element_fns.push_back(
560 _sub_elements[i]->template dof_transformation_right_fn<U>(ttype));
561 }
562
563 return [this, sub_element_fns](std::span<U> data,
564 std::span<const std::uint32_t> cell_info,
565 std::int32_t cell, int block_size)
566 {
567 std::size_t offset = 0;
568 for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
569 {
570 sub_element_fns[e](data.subspan(offset, data.size() - offset),
571 cell_info, cell, block_size);
572 offset += _sub_elements[e]->space_dimension();
573 }
574 };
575 }
576 else if (!scalar_element)
577 {
578 // Blocked element
579 // The transformation from the left can be used here as blocked
580 // elements use xyzxyzxyz ordering, and so applying the DOF
581 // transformation from the right is equivalent to applying the DOF
582 // transformation from the left to data using xxxyyyzzz ordering
583 std::function<void(std::span<U>, std::span<const std::uint32_t>,
584 std::int32_t, int)>
585 sub_fn
586 = _sub_elements.front()->template dof_transformation_fn<U>(ttype);
587 return [this, sub_fn](std::span<U> data,
588 std::span<const std::uint32_t> cell_info,
589 std::int32_t cell, int data_block_size)
590 {
591 const int ebs = block_size();
592 const std::size_t dof_count = data.size() / data_block_size;
593 for (int block = 0; block < data_block_size; ++block)
594 {
595 sub_fn(data.subspan(block * dof_count, dof_count), cell_info, cell,
596 ebs);
597 }
598 };
599 }
600 }
601
602 switch (ttype)
603 {
605 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
606 std::int32_t cell, int n)
607 { Tt_inv_apply_right(data, cell_info[cell], n); };
609 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
610 std::int32_t cell, int n)
611 { Tt_apply_right(data, cell_info[cell], n); };
613 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
614 std::int32_t cell, int n)
615 { Tinv_apply_right(data, cell_info[cell], n); };
617 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
618 std::int32_t cell, int n)
619 { T_apply_right(data, cell_info[cell], n); };
620 default:
621 throw std::runtime_error("Unknown transformation type");
622 }
623 }
624
661 template <typename U>
662 void T_apply(std::span<U> data, std::uint32_t cell_permutation, int n) const
663 {
664 assert(_element);
665 _element->T_apply(data, n, cell_permutation);
666 }
667
677 template <typename U>
678 void Tt_inv_apply(std::span<U> data, std::uint32_t cell_permutation,
679 int n) const
680 {
681 assert(_element);
682 _element->Tt_inv_apply(data, n, cell_permutation);
683 }
684
694 template <typename U>
695 void Tt_apply(std::span<U> data, std::uint32_t cell_permutation, int n) const
696 {
697 assert(_element);
698 _element->Tt_apply(data, n, cell_permutation);
699 }
700
709 template <typename U>
710 void Tinv_apply(std::span<U> data, std::uint32_t cell_permutation,
711 int n) const
712 {
713 assert(_element);
714 _element->Tinv_apply(data, n, cell_permutation);
715 }
716
725 template <typename U>
726 void T_apply_right(std::span<U> data, std::uint32_t cell_permutation,
727 int n) const
728 {
729 assert(_element);
730 _element->T_apply_right(data, n, cell_permutation);
731 }
732
742 template <typename U>
743 void Tinv_apply_right(std::span<U> data, std::uint32_t cell_permutation,
744 int n) const
745 {
746 assert(_element);
747 _element->Tinv_apply_right(data, n, cell_permutation);
748 }
749
759 template <typename U>
760 void Tt_apply_right(std::span<U> data, std::uint32_t cell_permutation,
761 int n) const
762 {
763 assert(_element);
764 _element->Tt_apply_right(data, n, cell_permutation);
765 }
766
776 template <typename U>
777 void Tt_inv_apply_right(std::span<U> data, std::uint32_t cell_permutation,
778 int n) const
779 {
780 assert(_element);
781 _element->Tt_inv_apply_right(data, n, cell_permutation);
782 }
783
800 void permute(std::span<std::int32_t> doflist,
801 std::uint32_t cell_permutation) const;
802
819 void permute_inv(std::span<std::int32_t> doflist,
820 std::uint32_t cell_permutation) const;
821
836 std::function<void(std::span<std::int32_t>, std::uint32_t)>
837 dof_permutation_fn(bool inverse = false, bool scalar_element = false) const;
838
839private:
840 // Value shape. For blocked elements this is larger than
841 // _reference_value_shape. For non-blocked 'primal' elements it is
842 // equal to _reference_value_shape. For mixed elements, it is
843 // std::nullopt.
844 std::optional<std::vector<std::size_t>> _value_shape;
845
846 // Block size for BlockedElements. This gives the number of DOFs
847 // co-located at each dof 'point'.
848 int _bs;
849
850 // Element cell shape
851 mesh::CellType _cell_type;
852
853 // Element signature
854 std::string _signature;
855
856 // Dimension of the finite element space (accounting for any blocking)
857 int _space_dim;
858
859 // List of sub-elements (if any)
860 std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>
861 _sub_elements;
862
863 // Value space shape, e.g. {} for a scalar, {3, 3} for a tensor in 3D.
864 // For a mixed element it is std::nullopt.
865 std::optional<std::vector<std::size_t>> _reference_value_shape;
866
867 // Basix Element (nullptr for mixed elements)
868 std::unique_ptr<basix::FiniteElement<geometry_type>> _element;
869
870 // Indicate whether this element represents a symmetric 2-tensor
871 bool _symmetric;
872
873 // Indicate whether the element needs permutations or transformations
874 bool _needs_dof_permutations;
875 bool _needs_dof_transformations;
876
877 std::vector<std::vector<std::vector<int>>> _entity_dofs;
878 std::vector<std::vector<std::vector<int>>> _entity_closure_dofs;
879
880 // Quadrature points of a quadrature element (0 dimensional array for
881 // all elements except quadrature elements)
882 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _points;
883};
884
885} // namespace dolfinx::fem
Definition CoordinateElement.h:26
bool symmetric() const
Does the element represent a symmetric 2-tensor?
Definition FiniteElement.cpp:353
FiniteElement & operator=(const FiniteElement &element)=delete
Copy assignment.
const std::vector< std::shared_ptr< const FiniteElement< geometry_type > > > & sub_elements() const noexcept
Get subelements (if any)
Definition FiniteElement.cpp:396
void Tinv_apply(std::span< U > data, std::uint32_t cell_permutation, int n) const
Apply the inverse of the operator applied by T_apply().
Definition FiniteElement.h:710
std::span< const std::size_t > value_shape() const
Value shape of the finite element field.
Definition FiniteElement.cpp:309
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > create_interpolation_operator(const FiniteElement &from) const
Create a matrix that maps degrees of freedom from one element to this element (interpolation).
Definition FiniteElement.cpp:496
void T_apply(std::span< U > data, std::uint32_t cell_permutation, int n) const
Transform basis functions from the reference element ordering and orientation to the globally consist...
Definition FiniteElement.h:662
basix::maps::type map_type() const
Get the map type used by the element.
Definition FiniteElement.cpp:425
int space_dimension() const noexcept
Dimension of the finite element function space (the number of degrees-of-freedom for the element).
Definition FiniteElement.cpp:291
void tabulate(std::span< geometry_type > values, std::span< const geometry_type > X, std::array< std::size_t, 2 > shape, int order) const
Evaluate derivatives of the basis functions up to given order at points in the reference cell.
Definition FiniteElement.cpp:365
void Tt_inv_apply_right(std::span< U > data, std::uint32_t cell_permutation, int n) const
Right(post)-apply the transpose inverse of the operator applied by T_apply().
Definition FiniteElement.h:777
FiniteElement & operator=(FiniteElement &&element)=default
Move assignment.
bool needs_dof_transformations() const noexcept
Check if DOF transformations are needed for this element.
Definition FiniteElement.cpp:539
FiniteElement(const FiniteElement &element)=delete
Copy constructor.
void Tt_inv_apply(std::span< U > data, std::uint32_t cell_permutation, int n) const
Apply the inverse transpose of the operator applied by T_apply().
Definition FiniteElement.h:678
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > interpolation_points() const
Points on the reference cell at which an expression needs to be evaluated in order to interpolate the...
Definition FiniteElement.cpp:464
std::shared_ptr< const FiniteElement< geometry_type > > extract_sub_element(const std::vector< int > &component) const
Extract sub finite element for component.
Definition FiniteElement.cpp:403
void T_apply_right(std::span< U > data, std::uint32_t cell_permutation, int n) const
Right(post)-apply the operator applied by T_apply().
Definition FiniteElement.h:726
std::string signature() const noexcept
String identifying the finite element.
Definition FiniteElement.cpp:285
bool needs_dof_permutations() const noexcept
Check if DOF permutations are needed for this element.
Definition FiniteElement.cpp:545
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const noexcept
Local DOFs associated with each sub-entity of the cell.
Definition FiniteElement.cpp:340
bool is_mixed() const noexcept
Check if element is a mixed element.
Definition FiniteElement.cpp:389
int reference_value_size() const
Value size of the base (non-blocked) finite element field.
Definition FiniteElement.cpp:318
int value_size() const
Value size of the finite element field.
Definition FiniteElement.cpp:297
std::function< void(std::span< U >, std::span< const std::uint32_t >, std::int32_t, int)> dof_transformation_fn(doftransform ttype, bool scalar_element=false) const
Return a function that applies a DOF transformation operator to some data (see T_apply()).
Definition FiniteElement.h:435
bool operator!=(const FiniteElement &e) const
Check if two elements are not equivalent.
Definition FiniteElement.cpp:273
~FiniteElement()=default
Destructor.
T geometry_type
Geometry type of the Mesh that the FunctionSpace is defined on.
Definition FiniteElement.h:60
const basix::FiniteElement< geometry_type > & basix_element() const
Return underlying Basix element (if it exists).
Definition FiniteElement.cpp:413
std::function< void(std::span< std::int32_t >, std::uint32_t)> dof_permutation_fn(bool inverse=false, bool scalar_element=false) const
Return a function that applies a degree-of-freedom permutation to some data.
FiniteElement(FiniteElement &&element)=default
Move constructor.
int num_sub_elements() const noexcept
Number of sub elements (for a mixed or blocked element).
Definition FiniteElement.cpp:383
void Tt_apply(std::span< U > data, std::uint32_t cell_permutation, int n) const
Apply the transpose of the operator applied by T_apply().
Definition FiniteElement.h:695
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const noexcept
Local DOFs associated with the closure of each sub-entity of the cell.
Definition FiniteElement.cpp:347
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > interpolation_operator() const
Definition FiniteElement.cpp:483
int block_size() const noexcept
Block size of the finite element function space.
Definition FiniteElement.cpp:359
std::span< const std::size_t > reference_value_shape() const
Value shape of the base (non-blocked) finite element field.
Definition FiniteElement.cpp:330
void Tt_apply_right(std::span< U > data, std::uint32_t cell_permutation, int n) const
Right(post)-apply the transpose of the operator applied by T_apply().
Definition FiniteElement.h:760
FiniteElement(const basix::FiniteElement< geometry_type > &element, std::optional< std::vector< std::size_t > > value_shape=std::nullopt, bool symmetric=false)
Create a finite element from a Basix finite element.
Definition FiniteElement.cpp:111
std::function< void(std::span< U >, std::span< const std::uint32_t >, std::int32_t, int)> dof_transformation_right_fn(doftransform ttype, bool scalar_element=false) const
Return a function that applies DOF transformation to some transposed data (see T_apply_right()).
Definition FiniteElement.h:539
mesh::CellType cell_type() const noexcept
Cell shape that the element is defined on.
Definition FiniteElement.cpp:279
bool interpolation_ident() const noexcept
Definition FiniteElement.cpp:451
void Tinv_apply_right(std::span< U > data, std::uint32_t cell_permutation, int n) const
Right(post)-apply the inverse of the operator applied by T_apply().
Definition FiniteElement.h:743
bool operator==(const FiniteElement &e) const
Check if two elements are equivalent.
Definition FiniteElement.cpp:261
bool map_ident() const noexcept
Definition FiniteElement.cpp:437
Finite element method functionality.
Definition assemble_expression_impl.h:23
BasixElementData(U element, V bs, W symmetry) -> BasixElementData< typename std::remove_cvref< U >::type::scalar_type >
Type deduction helper.
doftransform
DOF transformation type.
Definition FiniteElement.h:26
@ transpose
Transpose.
Definition FiniteElement.h:28
@ inverse_transpose
Transpose inverse.
Definition FiniteElement.h:30
@ inverse
Inverse.
Definition FiniteElement.h:29
@ standard
Standard.
Definition FiniteElement.h:27
@ cell
Cell.
Definition Form.h:37
CellType
Cell type identifier.
Definition cell_types.h:22
Basix element holder.
Definition FiniteElement.h:37
std::optional< std::vector< std::size_t > > value_shape
Value shape. Can only be set for scalar element.
Definition FiniteElement.h:41
bool symmetry
Definition FiniteElement.h:42
std::reference_wrapper< const basix::FiniteElement< T > > element
Finite element.
Definition FiniteElement.h:39