1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
//! Traits.
use crate::types::DofTransformation;
use rlst::{Array, MutableArrayImpl, RlstScalar, ValueArrayImpl};
use std::fmt::Debug;
use std::hash::Hash;
/// A finite element.
pub trait FiniteElement {
/// The scalar type
type T: RlstScalar;
/// Cell type
///
/// The cell type is typically defined through [ReferenceCellType](crate::types::ReferenceCellType).
type CellType: Debug + PartialEq + Eq + Clone + Copy + Hash;
/// The reference cell type, eg one of `Point`, `Interval`, `Triangle`, etc.
fn cell_type(&self) -> Self::CellType;
/// The number of basis functions.
fn dim(&self) -> usize;
/// The shape of the values returned by functions in $\mathcal{V}$.
///
/// If the values are scalar an empty slice is returned.
fn value_shape(&self) -> &[usize];
/// The number of values returned.
///
/// If (for example) `value_shape` is `[3, 4]` then `value_size` is $3\times4 = 12$.
/// If `value_shape` returns an empty array (ie the shape functions are scalar) the
// convention is used that the product of the elements of an empty array is 1.
fn value_size(&self) -> usize;
/// Tabulate the values of the basis functions and their derivatives at a set of points
///
/// - `points` is a two-dimensional array where each column contains the coordinates of one point.
/// - `nderivs` is the desired number of derivatives (0 for no derivatives, 1 for all first order derivatives, etc.).
/// - `data` is the 4-dimensional output array. The first dimension is the total number of partial derivatives,
/// the second dimension is the number of evaluation points, the third dimension is the number of
/// basis functions of the element, and the fourth dimension is the value size of the basis function output.
/// For example, `data[3, 2, 1, 0]` returns the 0th value of the third partial derivative on the point with index 2 for the
/// basis function with index 1.
///
/// ## Remark
/// Let $d^{i + k} = dx^{i}dy^{j}$ be a derivative with respect to $x$, $y$ in two dimensions and
/// $d^{i + k + j} = dx^{i}dy^{j}dz^{k}$ be a derivative with respect to $x$, $y$, and $z$ in three dimensions.
/// Then the corresponding index $\ell$ in the first dimension of the `data` array is computed as follows.
///
/// - Triangle: $l = (i + j + 1) * (i + j) / 2 + j$
/// - Quadrilateral: $l = i * (n + 1) + j$
/// - Tetrahedron: $l = (i + j + k) * (i + j + k + 1) * (i + j + k + 2) / 6 + (j + k) * (j + k + 1) / 2 + k$
/// - Hexahedron $l = i * (n + 1) * (n + 1) + j * (n + 1) + k$.
///
/// For the quadrilaterals and hexahedra, the parameter $n$ denotes the Lagrange superdegree.
fn tabulate<
TGeo: RlstScalar,
Array2Impl: ValueArrayImpl<TGeo, 2>,
Array4MutImpl: MutableArrayImpl<Self::T, 4>,
>(
&self,
points: &Array<Array2Impl, 2>,
nderivs: usize,
data: &mut Array<Array4MutImpl, 4>,
);
/// Get the required shape for a tabulation array.
fn tabulate_array_shape(&self, nderivs: usize, npoints: usize) -> [usize; 4];
/// Return the dof indices that are associated with the subentity with index `entity_number` and dimension `entity_dim`.
///
/// - For `entity_dim = 0` this returns the degrees of freedom (dofs) associated with the corresponding point.
/// - For `entity_dim = 1` this returns the dofs associated with the corresponding edge.
/// - For `entity_dim = 2` this returns the dofs associated with the corresponding face.
///
/// Note that this does not return dofs on the boundary of an entity, that means that (for example)
/// for an edge the dofs associated with the two vertices at the boundary of the edge are not returned.
/// To return also the boundary dofs use [FiniteElement::entity_closure_dofs].
fn entity_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>;
/// The DOFs that are associated with a closure of a subentity of the reference cell.
///
/// This method is similar to [FiniteElement::entity_dofs], but it includes the dofs
/// associated with the boundary of an entity. For an edge (for example) it returns the dofs associated
/// with the vertices at the boundary of the edge (as well as the dofs associated with the edge itself).
fn entity_closure_dofs(&self, entity_dim: usize, entity_number: usize) -> Option<&[usize]>;
}
/// A finite element that is mapped from a reference cell.
pub trait MappedFiniteElement: FiniteElement {
/// Transformation type
///
/// The Transformation type specifies possible transformations of the dofs on the reference element.
/// In most cases these will be rotations and reflections as defined in [Transformation](crate::types::Transformation).
type TransformationType: Debug + PartialEq + Eq + Clone + Copy + Hash;
/// The smallest degree Lagrange space that contains all possible polynomials of the finite element's polynomial space.
///
/// Details on the definition of the degree of Lagrange spaces of finite elements are
/// given [here](https://defelement.org/ciarlet.html#The+degree+of+a+finite+element).
fn lagrange_superdegree(&self) -> usize;
/// Push function values forward to a physical cell.
///
/// For Lagrange elements, this is an identity map. For many other element types, the function values
/// on the reference cell differ from those on the physical cell: for example Nedlec elements use a covariant
/// Piola transform. This method implements the appropriate transformation for the element.
///
/// - `reference_values`: The values on the reference cell. The shape of this input is the same as the `data` input to the function
/// [[FiniteElement::tabulate]].
/// - `nderivs`: The number of derivatives.
/// - `jacobians` should have shape [geometry_dimension, entity_topology_dimension, npts] and use column-major ordering;
/// - `jacobian_determinants` should have shape \[npts\];
/// - `inverse_jacobians` should have shape [entity_topology_dimension, geometry_dimension, npts] and use column-major ordering;
/// - `physical_values`: The output array of the push operation. This shape of this array is the same as the `reference_values`
/// but with reference value size replaced by the physical value size
fn push_forward<
TGeo: RlstScalar,
Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
Array4Impl: ValueArrayImpl<Self::T, 4>,
Array4MutImpl: MutableArrayImpl<Self::T, 4>,
>(
&self,
reference_values: &Array<Array4Impl, 4>,
nderivs: usize,
jacobians: &Array<Array3GeoImpl, 3>,
jacobian_determinants: &[TGeo],
inverse_jacobians: &Array<Array3GeoImpl, 3>,
physical_values: &mut Array<Array4MutImpl, 4>,
);
/// Pull function values back to the reference cell.
///
/// This is the inverse operation to [MappedFiniteElement::push_forward].
fn pull_back<
TGeo: RlstScalar,
Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
Array4Impl: ValueArrayImpl<Self::T, 4>,
Array4MutImpl: MutableArrayImpl<Self::T, 4>,
>(
&self,
physical_values: &Array<Array4Impl, 4>,
nderivs: usize,
jacobians: &Array<Array3GeoImpl, 3>,
jacobian_determinants: &[TGeo],
inverse_jacobians: &Array<Array3GeoImpl, 3>,
reference_values: &mut Array<Array4MutImpl, 4>,
);
/// The value shape on a physical cell
fn physical_value_shape(&self, gdim: usize) -> Vec<usize>;
/// The value size on a physical cell
fn physical_value_size(&self, gdim: usize) -> usize {
let mut vs = 1;
for i in self.physical_value_shape(gdim) {
vs *= i;
}
vs
}
/// The DOF transformation for a sub-entity.
fn dof_transformation(
&self,
entity: Self::CellType,
transformation: Self::TransformationType,
) -> Option<&DofTransformation<Self::T>>;
/// Apply permutation parts of DOF transformations.
fn apply_dof_permutations<T>(&self, data: &mut [T], cell_orientation: i32);
/// Apply non-permutation parts of DOF transformations.
fn apply_dof_transformations(&self, data: &mut [Self::T], cell_orientation: i32);
/// Apply DOF transformations.
fn apply_dof_permutations_and_transformations(
&self,
data: &mut [Self::T],
cell_orientation: i32,
) {
self.apply_dof_permutations(data, cell_orientation);
self.apply_dof_transformations(data, cell_orientation);
}
}
/// A family of finite elements defined on various cell types, with appropriate continuity
/// between different cell types.
pub trait ElementFamily {
/// The scalar type
type T: RlstScalar;
/// Cell type
type CellType: Debug + PartialEq + Eq + Clone + Copy + Hash;
/// The finite element type
type FiniteElement: FiniteElement<T = Self::T, CellType = Self::CellType> + 'static;
/// Create an element for the given cell type.
fn element(&self, cell_type: Self::CellType) -> Self::FiniteElement;
}
/// A map from the reference cell to physical cells.
pub trait Map {
/// Push values from the reference cell to physical cells.
fn push_forward<
T: RlstScalar,
TGeo: RlstScalar,
Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
Array4Impl: ValueArrayImpl<T, 4>,
Array4MutImpl: MutableArrayImpl<T, 4>,
>(
&self,
reference_values: &Array<Array4Impl, 4>,
nderivs: usize,
jacobians: &Array<Array3GeoImpl, 3>,
jacobian_determinants: &[TGeo],
inverse_jacobians: &Array<Array3GeoImpl, 3>,
physical_values: &mut Array<Array4MutImpl, 4>,
);
/// Pull values back to the reference cell.
fn pull_back<
T: RlstScalar,
TGeo: RlstScalar,
Array3GeoImpl: ValueArrayImpl<TGeo, 3>,
Array4Impl: ValueArrayImpl<T, 4>,
Array4MutImpl: MutableArrayImpl<T, 4>,
>(
&self,
physical_values: &Array<Array4Impl, 4>,
nderivs: usize,
jacobians: &Array<Array3GeoImpl, 3>,
jacobian_determinants: &[TGeo],
inverse_jacobians: &Array<Array3GeoImpl, 3>,
reference_values: &mut Array<Array4MutImpl, 4>,
);
/// The value shape on a physical cell.
fn physical_value_shape(&self, gdim: usize) -> Vec<usize>;
}