Skip to content

Commit

Permalink
add macro for derivative test
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Sep 2, 2024
1 parent 3481c88 commit d909add
Showing 1 changed file with 156 additions and 114 deletions.
270 changes: 156 additions & 114 deletions src/polynomials.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usi
#[cfg(test)]
mod test {
use super::*;
use crate::types::{Array2D, ReferenceCellType};
use crate::{quadrature::make_gauss_jacobi_quadrature, traits::QuadratureRule};
use approx::*;
use paste::paste;
Expand Down Expand Up @@ -100,130 +101,171 @@ mod test {
test_orthogonal!(Tetrahedron, 5);
test_orthogonal!(Tetrahedron, 6);

#[test]
fn test_legendre_interval_derivative() {
let degree = 6;

let epsilon = 1e-10;
let mut points = rlst_dynamic_array2!(f64, [1, 20]);
for i in 0..10 {
*points.get_mut([0, 2 * i]).unwrap() = i as f64 / 10.0;
*points.get_mut([0, 2 * i + 1]).unwrap() = points.get([0, 2 * i]).unwrap() + epsilon;
}

let mut data = rlst_dynamic_array3!(
f64,
legendre_shape(ReferenceCellType::Interval, &points, degree, 1,)
);
tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 1, &mut data);

for i in 0..degree + 1 {
for k in 0..points.shape()[1] / 2 {
assert_relative_eq!(
*data.get([1, i, 2 * k]).unwrap(),
(data.get([0, i, 2 * k + 1]).unwrap() - data.get([0, i, 2 * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
fn generate_points(cell: ReferenceCellType, epsilon: f64) -> Array2D<f64> {
match cell {
ReferenceCellType::Interval => {
let mut points = rlst_dynamic_array2!(f64, [1, 20]);
for i in 0..10 {
*points.get_mut([0, 2 * i]).unwrap() = i as f64 / 10.0;
*points.get_mut([0, 2 * i + 1]).unwrap() =
points.get([0, 2 * i]).unwrap() + epsilon;
}
points
}
}
}

#[test]
fn test_legendre_triangle_derivative() {
let degree = 6;

let epsilon = 1e-10;
let mut points = rlst_dynamic_array2!(f64, [2, 165]);
let mut index = 0;
for i in 0..10 {
for j in 0..10 - i {
*points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
*points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
*points.get_mut([0, 3 * index + 1]).unwrap() =
*points.get([0, 3 * index]).unwrap() + epsilon;
*points.get_mut([1, 3 * index + 1]).unwrap() = *points.get([1, 3 * index]).unwrap();
*points.get_mut([0, 3 * index + 2]).unwrap() = *points.get([0, 3 * index]).unwrap();
*points.get_mut([1, 3 * index + 2]).unwrap() =
*points.get([1, 3 * index]).unwrap() + epsilon;
index += 1;
ReferenceCellType::Triangle => {
let mut points = rlst_dynamic_array2!(f64, [2, 165]);
let mut index = 0;
for i in 0..10 {
for j in 0..10 - i {
*points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
*points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
*points.get_mut([0, 3 * index + 1]).unwrap() =
*points.get([0, 3 * index]).unwrap() + epsilon;
*points.get_mut([1, 3 * index + 1]).unwrap() =
*points.get([1, 3 * index]).unwrap();
*points.get_mut([0, 3 * index + 2]).unwrap() =
*points.get([0, 3 * index]).unwrap();
*points.get_mut([1, 3 * index + 2]).unwrap() =
*points.get([1, 3 * index]).unwrap() + epsilon;
index += 1;
}
}
points
}
}

let mut data = rlst_dynamic_array3!(
f64,
legendre_shape(ReferenceCellType::Triangle, &points, degree, 1,)
);
tabulate_legendre_polynomials(ReferenceCellType::Triangle, &points, degree, 1, &mut data);

for i in 0..degree + 1 {
for k in 0..points.shape()[1] / 3 {
assert_relative_eq!(
*data.get([1, i, 3 * k]).unwrap(),
(data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
assert_relative_eq!(
*data.get([2, i, 3 * k]).unwrap(),
(data.get([0, i, 3 * k + 2]).unwrap() - data.get([0, i, 3 * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
ReferenceCellType::Quadrilateral => {
let mut points = rlst_dynamic_array2!(f64, [2, 300]);
for i in 0..10 {
for j in 0..10 {
let index = 10 * i + j;
*points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
*points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
*points.get_mut([0, 3 * index + 1]).unwrap() =
*points.get([0, 3 * index]).unwrap() + epsilon;
*points.get_mut([1, 3 * index + 1]).unwrap() =
*points.get([1, 3 * index]).unwrap();
*points.get_mut([0, 3 * index + 2]).unwrap() =
*points.get([0, 3 * index]).unwrap();
*points.get_mut([1, 3 * index + 2]).unwrap() =
*points.get([1, 3 * index]).unwrap() + epsilon;
}
}
points
}
}
}

#[test]
fn test_legendre_quadrilateral_derivative() {
let degree = 6;

let epsilon = 1e-10;
let mut points = rlst_dynamic_array2!(f64, [2, 300]);
for i in 0..10 {
for j in 0..10 {
let index = 10 * i + j;
*points.get_mut([0, 3 * index]).unwrap() = i as f64 / 10.0;
*points.get_mut([1, 3 * index]).unwrap() = j as f64 / 10.0;
*points.get_mut([0, 3 * index + 1]).unwrap() =
*points.get([0, 3 * index]).unwrap() + epsilon;
*points.get_mut([1, 3 * index + 1]).unwrap() = *points.get([1, 3 * index]).unwrap();
*points.get_mut([0, 3 * index + 2]).unwrap() = *points.get([0, 3 * index]).unwrap();
*points.get_mut([1, 3 * index + 2]).unwrap() =
*points.get([1, 3 * index]).unwrap() + epsilon;
ReferenceCellType::Tetrahedron => {
let mut points = rlst_dynamic_array2!(f64, [2, 140]);
let mut index = 0;
for i in 0..5 {
for j in 0..5 - i {
for k in 0..5 - i - j {
*points.get_mut([0, 4 * index]).unwrap() = i as f64 / 5.0;
*points.get_mut([1, 4 * index]).unwrap() = j as f64 / 5.0;
*points.get_mut([2, 4 * index]).unwrap() = k as f64 / 5.0;
*points.get_mut([0, 4 * index + 1]).unwrap() =
*points.get([0, 4 * index]).unwrap() + epsilon;
*points.get_mut([1, 4 * index + 1]).unwrap() =
*points.get([1, 4 * index]).unwrap();
*points.get_mut([2, 4 * index + 1]).unwrap() =
*points.get([2, 4 * index]).unwrap();
*points.get_mut([0, 4 * index + 2]).unwrap() =
*points.get([0, 4 * index]).unwrap();
*points.get_mut([1, 4 * index + 2]).unwrap() =
*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([0, 4 * index]).unwrap();
*points.get_mut([1, 4 * index + 2]).unwrap() =
*points.get([1, 4 * index]).unwrap();
*points.get_mut([2, 4 * index + 2]).unwrap() =
*points.get([2, 4 * index]).unwrap() + epsilon;
index += 1;
}
}
}
points
}
ReferenceCellType::Hexahedron => {
let mut points = rlst_dynamic_array2!(f64, [2, 5000]);
for i in 0..5 {
for j in 0..5 {
for k in 0..5 {
let index = 100 * i + 10 * j + k;
*points.get_mut([0, 4 * index]).unwrap() = i as f64 / 5.0;
*points.get_mut([1, 4 * index]).unwrap() = j as f64 / 5.0;
*points.get_mut([2, 4 * index]).unwrap() = k as f64 / 5.0;
*points.get_mut([0, 4 * index + 1]).unwrap() =
*points.get([0, 4 * index]).unwrap() + epsilon;
*points.get_mut([1, 4 * index + 1]).unwrap() =
*points.get([1, 4 * index]).unwrap();
*points.get_mut([2, 4 * index + 1]).unwrap() =
*points.get([2, 4 * index]).unwrap();
*points.get_mut([0, 4 * index + 2]).unwrap() =
*points.get([0, 4 * index]).unwrap();
*points.get_mut([1, 4 * index + 2]).unwrap() =
*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([0, 4 * index]).unwrap();
*points.get_mut([1, 4 * index + 2]).unwrap() =
*points.get([1, 4 * index]).unwrap();
*points.get_mut([2, 4 * index + 2]).unwrap() =
*points.get([2, 4 * index]).unwrap() + epsilon;
}
}
}
points
}
_ => {
panic!("Unsupported cell type: {cell:?}");
}
}
}
macro_rules! test_orthogonal {
($cell:ident, $degree:expr) => {
paste! {
#[test]
fn [<test_legendre_derivatives_ $cell:lower _ $degree>]() {
let dim = reference_cell::dim(ReferenceCellType::[<$cell>]);
let epsilon = 1e-10;
let points = generate_points(ReferenceCellType::[<$cell>], epsilon);

let mut data = rlst_dynamic_array3!(
f64,
legendre_shape(ReferenceCellType::Quadrilateral, &points, degree, 1,)
);
tabulate_legendre_polynomials(
ReferenceCellType::Quadrilateral,
&points,
degree,
1,
&mut data,
);
let mut data = rlst_dynamic_array3!(
f64,
legendre_shape(ReferenceCellType::[<$cell>], &points, [<$degree>], 1,)
);
tabulate_legendre_polynomials(ReferenceCellType::[<$cell>], &points, [<$degree>], 1, &mut data);

for i in 0..degree + 1 {
for k in 0..points.shape()[1] / 3 {
assert_relative_eq!(
*data.get([1, i, 3 * k]).unwrap(),
(data.get([0, i, 3 * k + 1]).unwrap() - data.get([0, i, 3 * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
assert_relative_eq!(
*data.get([2, i, 3 * k]).unwrap(),
(data.get([0, i, 3 * k + 2]).unwrap() - data.get([0, i, 3 * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
for i in 0..[<$degree>] + 1 {
for k in 0..points.shape()[1] / (dim + 1) {
for d in 1..dim + 1 {
assert_relative_eq!(
*data.get([d, i, (dim + 1) * k]).unwrap(),
(data.get([0, i, (dim + 1) * k + d]).unwrap() - data.get([0, i, (dim + 1) * k]).unwrap())
/ epsilon,
epsilon = 1e-4
);
}
}
}
}
}
}
};
}

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]
fn test_legendre_interval_against_known_polynomials() {
let degree = 3;
Expand Down

0 comments on commit d909add

Please sign in to comment.