From 9d4b58574625e3a9f5a8756f515f00ca0625963f Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 3 Sep 2024 08:15:18 +0100 Subject: [PATCH] correct derivatives in tetrahedron polynomials --- src/polynomials.rs | 54 ++++++++++++++++++++++--------------- src/polynomials/legendre.rs | 31 ++++++++++++++------- 2 files changed, 54 insertions(+), 31 deletions(-) diff --git a/src/polynomials.rs b/src/polynomials.rs index 60246bb..2c641bc 100644 --- a/src/polynomials.rs +++ b/src/polynomials.rs @@ -152,7 +152,7 @@ mod test { points } ReferenceCellType::Tetrahedron => { - let mut points = rlst_dynamic_array2!(f64, [2, 140]); + let mut points = rlst_dynamic_array2!(f64, [3, 140]); let mut index = 0; for i in 0..5 { for j in 0..5 - i { @@ -172,11 +172,11 @@ mod test { *points.get([1, 4 * index]).unwrap() + epsilon; *points.get_mut([2, 4 * index + 2]).unwrap() = *points.get([2, 4 * index]).unwrap(); - *points.get_mut([0, 4 * index + 2]).unwrap() = + *points.get_mut([0, 4 * index + 3]).unwrap() = *points.get([0, 4 * index]).unwrap(); - *points.get_mut([1, 4 * index + 2]).unwrap() = + *points.get_mut([1, 4 * index + 3]).unwrap() = *points.get([1, 4 * index]).unwrap(); - *points.get_mut([2, 4 * index + 2]).unwrap() = + *points.get_mut([2, 4 * index + 3]).unwrap() = *points.get([2, 4 * index]).unwrap() + epsilon; index += 1; } @@ -185,7 +185,7 @@ mod test { points } ReferenceCellType::Hexahedron => { - let mut points = rlst_dynamic_array2!(f64, [2, 5000]); + let mut points = rlst_dynamic_array2!(f64, [3, 5000]); for i in 0..5 { for j in 0..5 { for k in 0..5 { @@ -205,11 +205,11 @@ mod test { *points.get([1, 4 * index]).unwrap() + epsilon; *points.get_mut([2, 4 * index + 2]).unwrap() = *points.get([2, 4 * index]).unwrap(); - *points.get_mut([0, 4 * index + 2]).unwrap() = + *points.get_mut([0, 4 * index + 3]).unwrap() = *points.get([0, 4 * index]).unwrap(); - *points.get_mut([1, 4 * index + 2]).unwrap() = + *points.get_mut([1, 4 * index + 3]).unwrap() = *points.get([1, 4 * index]).unwrap(); - *points.get_mut([2, 4 * index + 2]).unwrap() = + *points.get_mut([2, 4 * index + 3]).unwrap() = *points.get([2, 4 * index]).unwrap() + epsilon; } } @@ -221,7 +221,7 @@ mod test { } } } - macro_rules! test_orthogonal { + macro_rules! test_derivatives { ($cell:ident, $degree:expr) => { paste! { #[test] @@ -253,18 +253,30 @@ mod test { }; } - test_orthogonal!(Interval, 3); - test_orthogonal!(Interval, 4); - test_orthogonal!(Interval, 5); - test_orthogonal!(Interval, 6); - test_orthogonal!(Triangle, 3); - test_orthogonal!(Triangle, 4); - test_orthogonal!(Triangle, 5); - test_orthogonal!(Triangle, 6); - test_orthogonal!(Quadrilateral, 3); - test_orthogonal!(Quadrilateral, 4); - test_orthogonal!(Quadrilateral, 5); - test_orthogonal!(Quadrilateral, 6); + test_derivatives!(Interval, 1); + test_derivatives!(Interval, 2); + test_derivatives!(Interval, 3); + test_derivatives!(Interval, 4); + test_derivatives!(Interval, 5); + test_derivatives!(Interval, 6); + test_derivatives!(Triangle, 1); + test_derivatives!(Triangle, 2); + test_derivatives!(Triangle, 3); + test_derivatives!(Triangle, 4); + test_derivatives!(Triangle, 5); + test_derivatives!(Triangle, 6); + test_derivatives!(Quadrilateral, 1); + test_derivatives!(Quadrilateral, 2); + test_derivatives!(Quadrilateral, 3); + test_derivatives!(Quadrilateral, 4); + test_derivatives!(Quadrilateral, 5); + test_derivatives!(Quadrilateral, 6); + test_derivatives!(Tetrahedron, 1); + test_derivatives!(Tetrahedron, 2); + test_derivatives!(Tetrahedron, 3); + test_derivatives!(Tetrahedron, 4); + test_derivatives!(Tetrahedron, 5); + test_derivatives!(Tetrahedron, 6); #[test] fn test_legendre_interval_against_known_polynomials() { diff --git a/src/polynomials/legendre.rs b/src/polynomials/legendre.rs index 3cdfb19..49609a8 100644 --- a/src/polynomials/legendre.rs +++ b/src/polynomials/legendre.rs @@ -480,6 +480,9 @@ fn tabulate_tetrahedron< for p in 1..degree + 1 { let a = T::from(2 * p - 1).unwrap() / T::from(p).unwrap(); for i in 0..points.shape()[1] { + let d = *data + .get([tet_index(kx, ky, kz), tet_index(0, 0, p - 1), i]) + .unwrap(); *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() = (T::from(*points.get([0, i]).unwrap()).unwrap() @@ -488,9 +491,7 @@ fn tabulate_tetrahedron< + T::from(*points.get([2, i]).unwrap()).unwrap() - T::from(1.0).unwrap()) * a - * *data - .get([tet_index(kx, ky, kz), tet_index(0, 0, p - 1), i]) - .unwrap(); + * d; } if kx > 0 { for i in 0..points.shape()[1] { @@ -529,8 +530,10 @@ fn tabulate_tetrahedron< .unwrap(); *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) - .unwrap() -= (T::from(*points.get([1, i]).unwrap()).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap() + .unwrap() -= (T::from( + *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + ) + .unwrap() - T::from(1.0).unwrap()) .powi(2) * d @@ -544,8 +547,10 @@ fn tabulate_tetrahedron< *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() -= T::from(ky * 2).unwrap() - * (T::from(*points.get([1, i]).unwrap()).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap() + * (T::from( + *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + ) + .unwrap() - T::from(1.0).unwrap()) * d * (a - T::from(1.0).unwrap()); @@ -571,8 +576,10 @@ fn tabulate_tetrahedron< *data .get_mut([tet_index(kx, ky, kz), tet_index(0, 0, p), i]) .unwrap() -= T::from(kz * 2).unwrap() - * (T::from(*points.get([1, i]).unwrap()).unwrap() - + T::from(*points.get([2, i]).unwrap()).unwrap() + * (T::from( + *points.get([1, i]).unwrap() + *points.get([2, i]).unwrap(), + ) + .unwrap() - T::from(1.0).unwrap()) * d * (a - T::from(1.0).unwrap()); @@ -762,7 +769,11 @@ fn tabulate_tetrahedron< .get([tet_index(kx, ky, kz - 1), tet_index(r, q, p), i]) .unwrap(); *data - .get_mut([tet_index(kx, ky, kz), tet_index(r, q, p), i]) + .get_mut([ + tet_index(kx, ky, kz), + tet_index(r + 1, q, p), + i, + ]) .unwrap() += T::from(2 * kz).unwrap() * ar * d; } }