diff --git a/src/ciarlet.rs b/src/ciarlet.rs index e9beea4..ae6bb87 100644 --- a/src/ciarlet.rs +++ b/src/ciarlet.rs @@ -968,6 +968,20 @@ mod test { check_dofs(e); } + #[test] + fn test_nedelec_2_triangle() { + let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); + assert_eq!(e.value_size(), 2); + check_dofs(e); + } + + #[test] + fn test_nedelec_3_triangle() { + let e = nedelec::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); + assert_eq!(e.value_size(), 2); + check_dofs(e); + } + #[test] fn test_nedelec_1_quadrilateral() { let e = nedelec::create(ReferenceCellType::Quadrilateral, 1, Continuity::Standard); @@ -1003,20 +1017,6 @@ mod test { check_dofs(e); } - #[test] - fn test_nedelec_2_triangle() { - let e = nedelec::create::(ReferenceCellType::Triangle, 2, Continuity::Standard); - assert_eq!(e.value_size(), 2); - check_dofs(e); - } - - #[test] - fn test_nedelec_3_triangle() { - let e = nedelec::create::(ReferenceCellType::Triangle, 3, Continuity::Standard); - assert_eq!(e.value_size(), 2); - check_dofs(e); - } - #[test] fn test_nedelec_2_quadrilateral() { let e = nedelec::create::(ReferenceCellType::Quadrilateral, 2, Continuity::Standard); @@ -1052,6 +1052,67 @@ mod test { check_dofs(e); } + #[test] + fn test_nedelec_1_hexahedron() { + let e = nedelec::create(ReferenceCellType::Hexahedron, 1, Continuity::Standard); + assert_eq!(e.value_size(), 3); + let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6)); + let mut points = rlst_dynamic_array2!(f64, [3, 6]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 1.0; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 0.8; + *points.get_mut([0, 2]).unwrap() = 0.0; + *points.get_mut([1, 2]).unwrap() = 1.0; + *points.get_mut([2, 2]).unwrap() = 1.0; + *points.get_mut([0, 3]).unwrap() = 0.5; + *points.get_mut([1, 3]).unwrap() = 0.0; + *points.get_mut([2, 3]).unwrap() = 0.1; + *points.get_mut([0, 4]).unwrap() = 1.0; + *points.get_mut([1, 4]).unwrap() = 0.5; + *points.get_mut([2, 4]).unwrap() = 0.5; + *points.get_mut([0, 5]).unwrap() = 0.5; + *points.get_mut([1, 5]).unwrap() = 0.5; + *points.get_mut([2, 5]).unwrap() = 1.0; + e.tabulate(&points, 0, &mut data); + + for pt in 0..6 { + let x = *points.get([0, pt]).unwrap(); + let y = *points.get([1, pt]).unwrap(); + let z = *points.get([2, pt]).unwrap(); + for (i, basis_f) in [ + [(1.0 - y) * (1.0 - z), 0.0, 0.0], + [0.0, (1.0 - x) * (1.0 - z), 0.0], + [0.0, 0.0, (1.0 - x) * (1.0 - y)], + [0.0, x * (1.0 - z), 0.0], + [0.0, 0.0, x * (1.0 - y)], + [y * (1.0 - z), 0.0, 0.0], + [0.0, 0.0, (1.0 - x) * y], + [0.0, 0.0, x * y], + [(1.0 - y) * z, 0.0, 0.0], + [0.0, (1.0 - x) * z, 0.0], + [0.0, x * z, 0.0], + [y * z, 0.0, 0.0], + ].iter().enumerate() + { + for (d, value) in basis_f.iter().enumerate() { + println!("[{i},{d}] {} =?= {}", *data.get([0, pt, i, d]).unwrap(), value); + assert_relative_eq!(*data.get([0, pt, i, d]).unwrap(), value, epsilon = 1e-14); + } + } + } + check_dofs(e); + } + + #[test] + fn test_nedelec_2_hexahedron() { + let e = nedelec::create::(ReferenceCellType::Hexahedron, 2, Continuity::Standard); + assert_eq!(e.value_size(), 3); + check_dofs(e); + } + macro_rules! test_entity_closure_dofs_lagrange { ($cell:ident, $degree:expr) => { paste! { diff --git a/src/ciarlet/nedelec.rs b/src/ciarlet/nedelec.rs index f1d268a..2b52d2e 100644 --- a/src/ciarlet/nedelec.rs +++ b/src/ciarlet/nedelec.rs @@ -401,21 +401,21 @@ fn create_tp, T: RlstScalar + Matr for k in 0..pdim_edge { *wcoeffs .get_mut([ - 3 * (i * pdim_edge.pow(2) + j * pdim_edge + k), + i * pdim_edge.pow(2) + j * pdim_edge + k, 0, i * pdim_edge.pow(2) + j * pdim_edge + k, ]) .unwrap() = T::from(1.0).unwrap(); *wcoeffs .get_mut([ - 3 * (i * pdim_edge.pow(2) + j * pdim_edge + k) + 1, + pdim_edge.pow(2) * pdim_edge_minus1 + k * pdim_edge * pdim_edge_minus1 + i * pdim_edge + j, 1, k * pdim_edge.pow(2) + i * pdim_edge + j, ]) .unwrap() = T::from(1.0).unwrap(); *wcoeffs .get_mut([ - 3 * (i * pdim_edge.pow(2) + j * pdim_edge + k) + 2, + pdim_edge.pow(2) * pdim_edge_minus1 * 2 + j * pdim_edge * pdim_edge_minus1 + k * pdim_edge_minus1 + i, 2, j * pdim_edge.pow(2) + k * pdim_edge + i, ]) @@ -542,25 +542,25 @@ fn create_tp, T: RlstScalar + Matr m[3].push(rlst_dynamic_array3!(T, [0, tdim, 0])) } else { let interior_q = - gauss_jacobi_rule(ReferenceCellType::Quadrilateral, 2 * degree - 1).unwrap(); + gauss_jacobi_rule(ReferenceCellType::Hexahedron, 2 * degree - 1).unwrap(); let interior_pts_t = interior_q .points .iter() .map(|i| TReal::from(*i).unwrap()) .collect::>(); - let interior_pts = rlst_array_from_slice2!(&interior_pts_t, [2, interior_q.npoints]); + let interior_pts = rlst_array_from_slice2!(&interior_pts_t, [3, interior_q.npoints]); let mut interior_phi = rlst_dynamic_array3![ T, legendre_shape( - ReferenceCellType::Quadrilateral, + ReferenceCellType::Hexahedron, &interior_pts, degree - 1, 0 ) ]; tabulate_legendre_polynomials( - ReferenceCellType::Quadrilateral, + ReferenceCellType::Hexahedron, &interior_pts, degree - 1, 0, @@ -590,7 +590,7 @@ fn create_tp, T: RlstScalar + Matr mat[[index, 0, w_i]] = T::from(*wt).unwrap() * interior_phi[[ 0, - k * pdim_edge_minus1 * pdim_edge_minus1 + k * pdim_edge_minus1.pow(2) + j * pdim_edge_minus1 + i, w_i, @@ -598,7 +598,7 @@ fn create_tp, T: RlstScalar + Matr mat[[index + 1, 1, w_i]] = T::from(*wt).unwrap() * interior_phi[[ 0, - i * pdim_edge_minus1 * pdim_edge_minus1 + i * pdim_edge_minus1.pow(2) + k * pdim_edge_minus1 + j, w_i, @@ -606,7 +606,7 @@ fn create_tp, T: RlstScalar + Matr mat[[index + 2, 2, w_i]] = T::from(*wt).unwrap() * interior_phi[[ 0, - j * pdim_edge_minus1 * pdim_edge_minus1 + j * pdim_edge_minus1.pow(2) + i * pdim_edge_minus1 + k, w_i, @@ -620,6 +620,21 @@ fn create_tp, T: RlstScalar + Matr } } + use rlst::RawAccess; + + for i in 0..12 { + println!("{:?}", x[1][i].data()); + println!("{:?}", m[1][i].data()); + } + println!(); + for i in 0..6 { + println!("{:?}", x[2][i].data()); + println!("{:?}", m[2][i].data()); + } + println!(); + println!("{:?}", x[3][0].data()); + println!("{:?}", m[3][0].data()); + CiarletElement::create( "Nedelec (first kind)".to_string(), cell_type,