Skip to content

Commit

Permalink
working on nedelec on hexes
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Nov 26, 2024
1 parent b688c4b commit f5acbc5
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 24 deletions.
89 changes: 75 additions & 14 deletions src/ciarlet.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -968,6 +968,20 @@ mod test {
check_dofs(e);
}

#[test]
fn test_nedelec_2_triangle() {
let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
assert_eq!(e.value_size(), 2);
check_dofs(e);
}

#[test]
fn test_nedelec_3_triangle() {
let e = nedelec::create::<f64>(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);
Expand Down Expand Up @@ -1003,20 +1017,6 @@ mod test {
check_dofs(e);
}

#[test]
fn test_nedelec_2_triangle() {
let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 2, Continuity::Standard);
assert_eq!(e.value_size(), 2);
check_dofs(e);
}

#[test]
fn test_nedelec_3_triangle() {
let e = nedelec::create::<f64>(ReferenceCellType::Triangle, 3, Continuity::Standard);
assert_eq!(e.value_size(), 2);
check_dofs(e);
}

#[test]
fn test_nedelec_2_quadrilateral() {
let e = nedelec::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
Expand Down Expand Up @@ -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],

Check warning on line 1095 in src/ciarlet.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet.rs
[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::<f64>(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! {
Expand Down
35 changes: 25 additions & 10 deletions src/ciarlet/nedelec.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,21 +401,21 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + 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();

Check warning on line 408 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
*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();

Check warning on line 415 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
*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,
])
Expand Down Expand Up @@ -542,25 +542,25 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + 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::<Vec<_>>();
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]);

Check warning on line 552 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
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,
Expand Down Expand Up @@ -590,23 +590,23 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + Matr
mat[[index, 0, w_i]] = T::from(*wt).unwrap()

Check warning on line 590 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
* interior_phi[[
0,
k * pdim_edge_minus1 * pdim_edge_minus1
k * pdim_edge_minus1.pow(2)
+ j * pdim_edge_minus1
+ i,
w_i,
]];
mat[[index + 1, 1, w_i]] = T::from(*wt).unwrap()
* interior_phi[[

Check warning on line 599 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
0,
i * pdim_edge_minus1 * pdim_edge_minus1
i * pdim_edge_minus1.pow(2)
+ k * pdim_edge_minus1
+ j,
w_i,
]];
mat[[index + 2, 2, w_i]] = T::from(*wt).unwrap()
* interior_phi[[

Check warning on line 607 in src/ciarlet/nedelec.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,serde,strict") 

Diff in /home/runner/work/ndelement/ndelement/src/ciarlet/nedelec.rs
0,
j * pdim_edge_minus1 * pdim_edge_minus1
j * pdim_edge_minus1.pow(2)
+ i * pdim_edge_minus1
+ k,
w_i,
Expand All @@ -620,6 +620,21 @@ fn create_tp<TReal: RlstScalar<Real = TReal>, T: RlstScalar<Real = TReal> + 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,
Expand Down

0 comments on commit f5acbc5

Please sign in to comment.