Skip to content

Commit

Permalink
correct derivatives in tetrahedron polynomials
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Sep 3, 2024
1 parent d909add commit 9d4b585
Show file tree
Hide file tree
Showing 2 changed files with 54 additions and 31 deletions.
54 changes: 33 additions & 21 deletions src/polynomials.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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;
}
Expand All @@ -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 {
Expand All @@ -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;
}
}
Expand All @@ -221,7 +221,7 @@ mod test {
}
}
}
macro_rules! test_orthogonal {
macro_rules! test_derivatives {
($cell:ident, $degree:expr) => {
paste! {
#[test]
Expand Down Expand Up @@ -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() {
Expand Down
31 changes: 21 additions & 10 deletions src/polynomials/legendre.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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] {
Expand Down Expand Up @@ -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
Expand All @@ -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());
Expand All @@ -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());
Expand Down Expand Up @@ -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;
}
}
Expand Down

0 comments on commit 9d4b585

Please sign in to comment.