DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
FiniteElement.h
1// Copyright (C) 2020-2021 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 <span>
18#include <utility>
19#include <vector>
20
21namespace dolfinx::fem
22{
24enum class doftransform
25{
26 standard = 0,
27 transpose = 1,
28 inverse = 2,
30};
31
36template <std::floating_point T>
38{
39public:
41 using geometry_type = T;
42
48 std::size_t block_size, bool symmetric = false);
49
53 const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
54 elements);
55
62 FiniteElement(mesh::CellType cell_type, std::span<const geometry_type> points,
63 std::array<std::size_t, 2> pshape, std::size_t block_size,
64 bool symmetric = false);
65
67 FiniteElement(const FiniteElement& element) = delete;
68
70 FiniteElement(FiniteElement&& element) = default;
71
73 ~FiniteElement() = default;
74
76 FiniteElement& operator=(const FiniteElement& element) = delete;
77
79 FiniteElement& operator=(FiniteElement&& element) = default;
80
85 bool operator==(const FiniteElement& e) const;
86
91 bool operator!=(const FiniteElement& e) const;
92
99 std::string signature() const noexcept;
100
104 int space_dimension() const noexcept;
105
110 int block_size() const noexcept;
111
115 int reference_value_size() const;
116
118 std::span<const std::size_t> reference_value_shape() const noexcept;
119
121 const std::vector<std::vector<std::vector<int>>>&
122 entity_dofs() const noexcept;
123
125 const std::vector<std::vector<std::vector<int>>>&
126 entity_closure_dofs() const noexcept;
127
129 bool symmetric() const;
130
142 void tabulate(std::span<geometry_type> values,
143 std::span<const geometry_type> X,
144 std::array<std::size_t, 2> shape, int order) const;
145
155 std::pair<std::vector<geometry_type>, std::array<std::size_t, 4>>
156 tabulate(std::span<const geometry_type> X, std::array<std::size_t, 2> shape,
157 int order) const;
158
161 int num_sub_elements() const noexcept;
162
170 bool is_mixed() const noexcept;
171
173 const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>&
174 sub_elements() const noexcept;
175
177 std::shared_ptr<const FiniteElement<geometry_type>>
178 extract_sub_element(const std::vector<int>& component) const;
179
182 const basix::FiniteElement<geometry_type>& basix_element() const;
183
185 basix::maps::type map_type() const;
186
193 bool interpolation_ident() const noexcept;
194
199 bool map_ident() const noexcept;
200
211 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
212 interpolation_points() const;
213
223 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
225
238 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
240
254 bool needs_dof_transformations() const noexcept;
255
270 bool needs_dof_permutations() const noexcept;
271
311 template <typename U>
312 std::function<void(std::span<U>, std::span<const std::uint32_t>, std::int32_t,
313 int)>
314 dof_transformation_fn(doftransform ttype, bool scalar_element = false) const
315 {
317 {
318 // If no permutation needed, return function that does nothing
319 return [](std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)
320 {
321 // Do nothing
322 };
323 }
324
325 if (!_sub_elements.empty())
326 {
327 if (_is_mixed)
328 {
329 // Mixed element
330 std::vector<std::function<void(
331 std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
332 sub_element_fns;
333 std::vector<int> dims;
334 for (std::size_t i = 0; i < _sub_elements.size(); ++i)
335 {
336 sub_element_fns.push_back(
337 _sub_elements[i]->template dof_transformation_fn<U>(ttype));
338 dims.push_back(_sub_elements[i]->space_dimension());
339 }
340
341 return [dims, sub_element_fns](std::span<U> data,
342 std::span<const std::uint32_t> cell_info,
343 std::int32_t cell, int block_size)
344 {
345 std::size_t offset = 0;
346 for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
347 {
348 const std::size_t width = dims[e] * block_size;
349 sub_element_fns[e](data.subspan(offset, width), cell_info, cell,
350 block_size);
351 offset += width;
352 }
353 };
354 }
355 else if (!scalar_element)
356 {
357 // Blocked element
358 std::function<void(std::span<U>, std::span<const std::uint32_t>,
359 std::int32_t, int)>
360 sub_fn = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
361 const int ebs = _bs;
362 return [ebs, sub_fn](std::span<U> data,
363 std::span<const std::uint32_t> cell_info,
364 std::int32_t cell, int data_block_size)
365 { sub_fn(data, cell_info, cell, ebs * data_block_size); };
366 }
367 }
368
369 switch (ttype)
370 {
372 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
373 std::int32_t cell, int block_size)
374 { Tt_inv_apply(data, cell_info[cell], block_size); };
376 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
377 std::int32_t cell, int block_size)
378 { Tt_apply(data, cell_info[cell], block_size); };
380 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
381 std::int32_t cell, int block_size)
382 { Tinv_apply(data, cell_info[cell], block_size); };
384 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
385 std::int32_t cell, int block_size)
386 { T_apply(data, cell_info[cell], block_size); };
387 default:
388 throw std::runtime_error("Unknown transformation type");
389 }
390 }
391
415 template <typename U>
416 std::function<void(std::span<U>, std::span<const std::uint32_t>, std::int32_t,
417 int)>
419 bool scalar_element = false) const
420 {
421 if (!needs_dof_transformations())
422 {
423 // If no permutation needed, return function that does nothing
424 return [](std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)
425 {
426 // Do nothing
427 };
428 }
429 else if (_sub_elements.size() != 0)
430 {
431 if (_is_mixed)
432 {
433 // Mixed element
434 std::vector<std::function<void(
435 std::span<U>, std::span<const std::uint32_t>, std::int32_t, int)>>
436 sub_element_fns;
437 for (std::size_t i = 0; i < _sub_elements.size(); ++i)
438 {
439 sub_element_fns.push_back(
440 _sub_elements[i]->template dof_transformation_right_fn<U>(ttype));
441 }
442
443 return [this, sub_element_fns](std::span<U> data,
444 std::span<const std::uint32_t> cell_info,
445 std::int32_t cell, int block_size)
446 {
447 std::size_t offset = 0;
448 for (std::size_t e = 0; e < sub_element_fns.size(); ++e)
449 {
450 sub_element_fns[e](data.subspan(offset, data.size() - offset),
451 cell_info, cell, block_size);
452 offset += _sub_elements[e]->space_dimension();
453 }
454 };
455 }
456 else if (!scalar_element)
457 {
458 // Blocked element
459 // The transformation from the left can be used here as blocked
460 // elements use xyzxyzxyz ordering, and so applying the DOF
461 // transformation from the right is equivalent to applying the DOF
462 // transformation from the left to data using xxxyyyzzz ordering
463 std::function<void(std::span<U>, std::span<const std::uint32_t>,
464 std::int32_t, int)>
465 sub_fn = _sub_elements[0]->template dof_transformation_fn<U>(ttype);
466 return [this, sub_fn](std::span<U> data,
467 std::span<const std::uint32_t> cell_info,
468 std::int32_t cell, int data_block_size)
469 {
470 const int ebs = block_size();
471 const std::size_t dof_count = data.size() / data_block_size;
472 for (int block = 0; block < data_block_size; ++block)
473 {
474 sub_fn(data.subspan(block * dof_count, dof_count), cell_info, cell,
475 ebs);
476 }
477 };
478 }
479 }
480
481 switch (ttype)
482 {
484 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
485 std::int32_t cell, int n)
486 { Tt_inv_apply_right(data, cell_info[cell], n); };
488 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
489 std::int32_t cell, int n)
490 { Tt_apply_right(data, cell_info[cell], n); };
492 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
493 std::int32_t cell, int n)
494 { Tinv_apply_right(data, cell_info[cell], n); };
496 return [this](std::span<U> data, std::span<const std::uint32_t> cell_info,
497 std::int32_t cell, int n)
498 { T_apply_right(data, cell_info[cell], n); };
499 default:
500 throw std::runtime_error("Unknown transformation type");
501 }
502 }
503
540 template <typename U>
541 void T_apply(std::span<U> data, std::uint32_t cell_permutation, int n) const
542 {
543 assert(_element);
544 _element->T_apply(data, n, cell_permutation);
545 }
546
556 template <typename U>
557 void Tt_inv_apply(std::span<U> data, std::uint32_t cell_permutation,
558 int n) const
559 {
560 assert(_element);
561 _element->Tt_inv_apply(data, n, cell_permutation);
562 }
563
573 template <typename U>
574 void Tt_apply(std::span<U> data, std::uint32_t cell_permutation, int n) const
575 {
576 assert(_element);
577 _element->Tt_apply(data, n, cell_permutation);
578 }
579
588 template <typename U>
589 void Tinv_apply(std::span<U> data, std::uint32_t cell_permutation,
590 int n) const
591 {
592 assert(_element);
593 _element->Tinv_apply(data, n, cell_permutation);
594 }
595
604 template <typename U>
605 void T_apply_right(std::span<U> data, std::uint32_t cell_permutation,
606 int n) const
607 {
608 assert(_element);
609 _element->T_apply_right(data, n, cell_permutation);
610 }
611
621 template <typename U>
622 void Tinv_apply_right(std::span<U> data, std::uint32_t cell_permutation,
623 int n) const
624 {
625 assert(_element);
626 _element->Tinv_apply_right(data, n, cell_permutation);
627 }
628
638 template <typename U>
639 void Tt_apply_right(std::span<U> data, std::uint32_t cell_permutation,
640 int n) const
641 {
642 assert(_element);
643 _element->Tt_apply_right(data, n, cell_permutation);
644 }
645
655 template <typename U>
656 void Tt_inv_apply_right(std::span<U> data, std::uint32_t cell_permutation,
657 int n) const
658 {
659 assert(_element);
660 _element->Tt_inv_apply_right(data, n, cell_permutation);
661 }
662
679 void permute(std::span<std::int32_t> doflist,
680 std::uint32_t cell_permutation) const;
681
698 void permute_inv(std::span<std::int32_t> doflist,
699 std::uint32_t cell_permutation) const;
700
715 std::function<void(std::span<std::int32_t>, std::uint32_t)>
716 dof_permutation_fn(bool inverse = false, bool scalar_element = false) const;
717
718private:
719 std::string _signature;
720
721 int _space_dim;
722
723 // List of sub-elements (if any)
724 std::vector<std::shared_ptr<const FiniteElement<geometry_type>>>
725 _sub_elements;
726
727 // Dimension of each value space
728 std::vector<std::size_t> _reference_value_shape;
729
730 // Block size for BlockedElements. This gives the number of DOFs
731 // co-located at each dof 'point'.
732 int _bs;
733
734 // Indicate whether this is a mixed element
735 bool _is_mixed;
736
737 // Basix Element (nullptr for mixed elements)
738 std::unique_ptr<basix::FiniteElement<geometry_type>> _element;
739
740 // Indicate whether this element represents a symmetric 2-tensor
741 bool _symmetric;
742
743 // Indicate whether the element needs permutations or transformations
744 bool _needs_dof_permutations;
745 bool _needs_dof_transformations;
746
747 std::vector<std::vector<std::vector<int>>> _entity_dofs;
748 std::vector<std::vector<std::vector<int>>> _entity_closure_dofs;
749
750 // Quadrature points of a quadrature element (0 dimensional array for
751 // all elements except quadrature elements)
752 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _points;
753};
754} // namespace dolfinx::fem
Definition utils.h:42
Model of a finite element.
Definition FiniteElement.h:38
bool symmetric() const
Does the element represent a symmetric 2-tensor?
Definition FiniteElement.cpp:250
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:300
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:589
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:390
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:541
basix::maps::type map_type() const
Get the map type used by the element.
Definition FiniteElement.cpp:329
int space_dimension() const noexcept
Definition FiniteElement.cpp:223
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:269
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:656
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:433
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:557
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:361
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:307
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:605
std::string signature() const noexcept
Definition FiniteElement.cpp:217
bool needs_dof_permutations() const noexcept
Check if DOF permutations are needed for this element.
Definition FiniteElement.cpp:439
const std::vector< std::vector< std::vector< int > > > & entity_dofs() const noexcept
The local DOFs associated with each subentity of the cell.
Definition FiniteElement.cpp:237
bool is_mixed() const noexcept
Check if element is a mixed element.
Definition FiniteElement.cpp:293
int reference_value_size() const
Definition FiniteElement.cpp:256
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:314
bool operator!=(const FiniteElement &e) const
Definition FiniteElement.cpp:211
~FiniteElement()=default
Destructor.
T geometry_type
Geometry type of the Mesh that the FunctionSpace is defined on.
Definition FiniteElement.h:41
const basix::FiniteElement< geometry_type > & basix_element() const
Return underlying Basix element (if it exists).
Definition FiniteElement.cpp:317
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(const basix::FiniteElement< geometry_type > &element, std::size_t block_size, bool symmetric=false)
Create a finite element from a Basix finite element.
Definition FiniteElement.cpp:97
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:287
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:574
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs() const noexcept
The local DOFs associated with the closure of each subentity of the cell.
Definition FiniteElement.cpp:244
std::span< const std::size_t > reference_value_shape() const noexcept
The reference value shape.
Definition FiniteElement.cpp:230
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > interpolation_operator() const
Definition FiniteElement.cpp:377
int block_size() const noexcept
Definition FiniteElement.cpp:263
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:639
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:418
bool interpolation_ident() const noexcept
Definition FiniteElement.cpp:351
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:622
bool operator==(const FiniteElement &e) const
Definition FiniteElement.cpp:199
bool map_ident() const noexcept
Definition FiniteElement.cpp:341
Finite element method functionality.
Definition assemble_matrix_impl.h:26
doftransform
DOF transformation type.
Definition FiniteElement.h:25
@ inverse_transpose
Transpose inverse.
CellType
Cell type identifier.
Definition cell_types.h:22