Skip to content

Commit 3e021f2

Browse files
mscroggsmscroggs
authored
Merge branch 'main' into mscroggs/system
2 parents df8555b + 0137405 commit 3e021f2

26 files changed

Lines changed: 1180 additions & 277 deletions

File tree

‎ndelement/Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ bempp-quadrature = { version = "0.1.0" }
2626
itertools = "0.14.*"
2727
mpi = { version = "0.8.0", optional = true }
2828
num = "0.4"
29-
rlst = "0.5"
29+
rlst = "0.6"
3030
serde = { version = "1", features = ["derive"], optional = true }
3131
strum = "0.27"
3232
strum_macros = "0.27"

‎ndelement/examples/lagrange_element.rs‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ fn main() {
2222
element.tabulate(&points, 0, &mut basis_values);
2323
println!(
2424
"The values of the basis functions at the point (1/3, 1/3) are: {:?}",
25-
basis_values.data()
25+
basis_values.data().unwrap()
2626
);
2727

2828
// Set point to [1, 0]
@@ -32,6 +32,6 @@ fn main() {
3232
element.tabulate(&points, 0, &mut basis_values);
3333
println!(
3434
"The values of the basis functions at the point (1, 0) are: {:?}",
35-
basis_values.data()
35+
basis_values.data().unwrap()
3636
);
3737
}

‎ndelement/python/test/test_element.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -273,17 +273,17 @@ def test_nedelec_1_triangle(continuity):
273273
def test_push_forward_pull_back(ftype, reference_values, physical_values):
274274
family = create_family(ftype, 1, continuity=Continuity.Standard)
275275
element = family.element(ReferenceCellType.Triangle)
276-
j = np.array([[[1.0, 2.0], [0.0, 0.0]], [[1.0, 0.0], [1.0, 3.0]]])
276+
j = np.array([[[1.0, 0.0], [1.0, 1.0]], [[2.0, 0.0], [0.0, 3.0]]])
277277
jdet = np.array([1.0, 6.0])
278278
jinv = np.array(
279279
[
280280
[
281-
[1.0, 0.5],
282-
[0.0, 0.0],
281+
[1.0, 0.0],
282+
[-1.0, 1.0],
283283
],
284284
[
285-
[-1.0, 0.0],
286-
[1.0, 1.0 / 3.0],
285+
[0.5, 0.0],
286+
[0.0, 1.0 / 3.0],
287287
],
288288
]
289289
)

‎ndelement/src/ciarlet.rs‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,7 @@ fn compute_derivative_count(nderivs: usize, cell_type: ReferenceCellType) -> usi
6565
}
6666

6767
/// A Ciarlet element
68-
pub struct CiarletElement<T: RlstScalar, M: Map, TGeo: RlstScalar = <T as RlstScalar>::Real> {
68+
pub struct CiarletElement<T: RlstScalar, M: Map, TGeo: RlstScalar> {
6969
family_name: String,
7070
cell_type: ReferenceCellType,
7171
degree: usize,
@@ -453,7 +453,7 @@ where
453453
let wts = &new_wts[edim][entity_id];
454454
let edofs = &entity_dofs[edim][entity_id];
455455
let endofs = edofs.len();
456-
let mut j = rlst_dynamic_array![TGeo, [npts, tdim, tdim]];
456+
let mut j = rlst_dynamic_array![TGeo, [tdim, tdim, npts]];
457457
for t_in in 0..tdim {
458458
for (t_out, (a, b)) in izip!(
459459
f(&vec![TGeo::zero(); tdim]),
@@ -466,7 +466,7 @@ where
466466
.enumerate()
467467
{
468468
for p in 0..npts {
469-
*j.get_mut([p, t_out, t_in]).unwrap() = b - a;
469+
*j.get_mut([t_out, t_in, p]).unwrap() = b - a;
470470
}
471471
}
472472
}
@@ -478,7 +478,7 @@ where
478478
};
479479
npts
480480
];
481-
let mut jinv = rlst_dynamic_array![TGeo, [npts, tdim, tdim]];
481+
let mut jinv = rlst_dynamic_array![TGeo, [tdim, tdim, npts]];
482482
for t_in in 0..tdim {
483483
for (t_out, (a, b)) in izip!(
484484
finv(&vec![TGeo::zero(); tdim], f),
@@ -494,7 +494,7 @@ where
494494
.enumerate()
495495
{
496496
for p in 0..npts {
497-
*jinv.get_mut([p, t_out, t_in]).unwrap() = TGeo::from(b - a).unwrap();
497+
*jinv.get_mut([t_out, t_in, p]).unwrap() = TGeo::from(b - a).unwrap();
498498
}
499499
}
500500
}

‎ndelement/src/map.rs‎

Lines changed: 23 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -85,8 +85,8 @@ impl Map for CovariantPiolaMap {
8585
inverse_jacobians: &Array<Array3GeoImpl, 3>,
8686
physical_values: &mut Array<Array4MutImpl, 4>,
8787
) {
88-
let tdim = inverse_jacobians.shape()[1];
89-
let gdim = inverse_jacobians.shape()[2];
88+
let tdim = inverse_jacobians.shape()[0];
89+
let gdim = inverse_jacobians.shape()[1];
9090
assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
9191
assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
9292
assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
@@ -101,7 +101,7 @@ impl Map for CovariantPiolaMap {
101101
unsafe {
102102
*physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim)
103103
.map(|td| {
104-
T::from(inverse_jacobians.get_value_unchecked([p, td, gd])).unwrap()
104+
T::from(inverse_jacobians.get_value_unchecked([td, gd, p])).unwrap()
105105
* reference_values.get_value_unchecked([0, p, b, td])
106106
})
107107
.sum::<T>();
@@ -125,8 +125,8 @@ impl Map for CovariantPiolaMap {
125125
_inverse_jacobians: &Array<Array3GeoImpl, 3>,
126126
reference_values: &mut Array<Array4MutImpl, 4>,
127127
) {
128-
let gdim = jacobians.shape()[1];
129-
let tdim = jacobians.shape()[2];
128+
let gdim = jacobians.shape()[0];
129+
let tdim = jacobians.shape()[1];
130130
assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
131131
assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
132132
assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
@@ -141,7 +141,7 @@ impl Map for CovariantPiolaMap {
141141
unsafe {
142142
*reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim)
143143
.map(|gd| {
144-
T::from(jacobians.get_value_unchecked([p, gd, td])).unwrap()
144+
T::from(jacobians.get_value_unchecked([gd, td, p])).unwrap()
145145
* physical_values.get_value_unchecked([0, p, b, gd])
146146
})
147147
.sum::<T>();
@@ -184,8 +184,9 @@ impl Map for ContravariantPiolaMap {
184184
_inverse_jacobians: &Array<Array3GeoImpl, 3>,
185185
physical_values: &mut Array<Array4MutImpl, 4>,
186186
) {
187-
let gdim = jacobians.shape()[1];
188-
let tdim = jacobians.shape()[2];
187+
let gdim = jacobians.shape()[0];
188+
let tdim = jacobians.shape()[1];
189+
189190
assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
190191
assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
191192
assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
@@ -200,7 +201,7 @@ impl Map for ContravariantPiolaMap {
200201
unsafe {
201202
*physical_values.get_unchecked_mut([0, p, b, gd]) = (0..tdim)
202203
.map(|td| {
203-
T::from(jacobians.get_value_unchecked([p, gd, td])).unwrap()
204+
T::from(jacobians.get_value_unchecked([gd, td, p])).unwrap()
204205
* reference_values.get_value_unchecked([0, p, b, td])
205206
})
206207
.sum::<T>()
@@ -225,8 +226,8 @@ impl Map for ContravariantPiolaMap {
225226
inverse_jacobians: &Array<Array3GeoImpl, 3>,
226227
reference_values: &mut Array<Array4MutImpl, 4>,
227228
) {
228-
let tdim = inverse_jacobians.shape()[1];
229-
let gdim = inverse_jacobians.shape()[2];
229+
let tdim = inverse_jacobians.shape()[0];
230+
let gdim = inverse_jacobians.shape()[1];
230231
assert_eq!(reference_values.shape()[0], physical_values.shape()[0]);
231232
assert_eq!(reference_values.shape()[1], physical_values.shape()[1]);
232233
assert_eq!(reference_values.shape()[2], physical_values.shape()[2]);
@@ -241,7 +242,7 @@ impl Map for ContravariantPiolaMap {
241242
unsafe {
242243
*reference_values.get_unchecked_mut([0, p, b, td]) = (0..gdim)
243244
.map(|gd| {
244-
T::from(inverse_jacobians.get_value_unchecked([p, td, gd])).unwrap()
245+
T::from(inverse_jacobians.get_value_unchecked([td, gd, p])).unwrap()
245246
* physical_values.get_value_unchecked([0, p, b, gd])
246247
})
247248
.sum::<T>()
@@ -268,22 +269,22 @@ mod test {
268269
jinv: &mut Array<Array3MutImpl, 3>,
269270
) {
270271
*j.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
271-
*j.get_mut([0, 0, 1]).unwrap() = T::from(1.0).unwrap();
272-
*j.get_mut([0, 1, 0]).unwrap() = T::from(0.0).unwrap();
273-
*j.get_mut([0, 1, 1]).unwrap() = T::from(1.0).unwrap();
274-
*j.get_mut([1, 0, 0]).unwrap() = T::from(2.0).unwrap();
272+
*j.get_mut([0, 1, 0]).unwrap() = T::from(1.0).unwrap();
273+
*j.get_mut([1, 0, 0]).unwrap() = T::from(0.0).unwrap();
274+
*j.get_mut([1, 1, 0]).unwrap() = T::from(1.0).unwrap();
275+
*j.get_mut([0, 0, 1]).unwrap() = T::from(2.0).unwrap();
276+
*j.get_mut([0, 1, 1]).unwrap() = T::from(0.0).unwrap();
275277
*j.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
276-
*j.get_mut([1, 1, 0]).unwrap() = T::from(0.0).unwrap();
277278
*j.get_mut([1, 1, 1]).unwrap() = T::from(3.0).unwrap();
278279
jdet[0] = T::from(1.0).unwrap();
279280
jdet[1] = T::from(6.0).unwrap();
280281
*jinv.get_mut([0, 0, 0]).unwrap() = T::from(1.0).unwrap();
281-
*jinv.get_mut([0, 0, 1]).unwrap() = T::from(-1.0).unwrap();
282-
*jinv.get_mut([0, 1, 0]).unwrap() = T::from(0.0).unwrap();
283-
*jinv.get_mut([0, 1, 1]).unwrap() = T::from(1.0).unwrap();
284-
*jinv.get_mut([1, 0, 0]).unwrap() = T::from(0.5).unwrap();
282+
*jinv.get_mut([0, 1, 0]).unwrap() = T::from(-1.0).unwrap();
283+
*jinv.get_mut([1, 0, 0]).unwrap() = T::from(0.0).unwrap();
284+
*jinv.get_mut([1, 1, 0]).unwrap() = T::from(1.0).unwrap();
285+
*jinv.get_mut([0, 0, 1]).unwrap() = T::from(0.5).unwrap();
286+
*jinv.get_mut([0, 1, 1]).unwrap() = T::from(0.0).unwrap();
285287
*jinv.get_mut([1, 0, 1]).unwrap() = T::from(0.0).unwrap();
286-
*jinv.get_mut([1, 1, 0]).unwrap() = T::from(0.0).unwrap();
287288
*jinv.get_mut([1, 1, 1]).unwrap() = T::from(1.0 / 3.0).unwrap();
288289
}
289290

‎ndelement/src/traits.rs‎

Lines changed: 5 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -107,21 +107,13 @@ pub trait MappedFiniteElement: FiniteElement {
107107
/// Piola transform. This method implements the appropriate transformation for the element.
108108
///
109109
/// - `reference_values`: The values on the reference cell. The shape of this input is the same as the `data` input to the function
110-
/// [[FiniteElement::tabulate].
110+
/// [[FiniteElement::tabulate]].
111111
/// - `nderivs`: The number of derivatives.
112-
/// - `jacobians:` A three-dimensional array of jacobians of the map from reference to physical cell.
113-
/// The first dimension is the reference point, the second dimension is the geometric dimension of the physical space, and
114-
/// the third dimension is the topological dimension of the reference element. For example,
115-
/// for the map of 5 points from the reference triangle to a physical surface triangle embedded in 3d space the dimension
116-
/// of `jacobians` is `[5, 3, 2]`.
117-
/// - `jacobian_determinants`: The determinant of the jacobian at each point. If the jacobian $J$ is not square, then the
118-
/// determinant is computed using $d=\sqrt{\det(J^TJ)}$.
119-
/// - `inverse_jacobians`: A three dimensional array that stores the inverse jacobian for the point with index j at position
120-
/// `inverse_jacobians[j, :, :]`. The first dimension of `inverse_jacobians` is the point index, the second dimension
121-
/// is the topological dimension, and the third dimension is the geometric dimension. If the Jacobian is rectangular then the
122-
/// inverse Jacobian is the pseudo-inverse of the Jacobian, ie the matrix $J^\dagger$ such that $J^\dagger J = I$.
112+
/// - `jacobians` should have shape [geometry_dimension, entity_topology_dimension, npts] and use column-major ordering;
113+
/// - `jacobian_determinants` should have shape \[npts\];
114+
/// - `inverse_jacobians` should have shape [entity_topology_dimension, geometry_dimension, npts] and use column-major ordering;
123115
/// - `physical_values`: The output array of the push operation. This shape of this array is the same as the `reference_values`
124-
/// input, with the [MappedFiniteElement::physical_value_size] used instead of the reference value size.
116+
/// but with reference value size replaced by the physical value size
125117
fn push_forward<
126118
TGeo: RlstScalar,
127119
Array3GeoImpl: ValueArrayImpl<TGeo, 3>,

‎ndfunctionspace/Cargo.toml‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,12 @@ ndelement = { path = "../ndelement" }
2626
ndgrid = { path = "../ndgrid" }
2727
itertools = "0.14.*"
2828
mpi = { version = "0.8.0", optional = true }
29-
rlst = { version = "0.5", optional = true }
29+
rlst = { version = "0.6", optional = true }
30+
31+
[dev-dependencies]
32+
quadraturerules = "0.9.0"
33+
approx = "0.5"
34+
rlst = "0.6"
3035

3136
[package.metadata.docs.rs]
3237
cargo-args = ["-Zunstable-options", "-Zrustdoc-scrape-examples"]

0 commit comments

Comments
 (0)