Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transpose point storage #10

Merged
merged 3 commits into from
Jul 25, 2024
Merged

Transpose point storage #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8afdc7d
start transposing point storage
mscroggs Jul 24, 2024
866c9a1
transpose points
mscroggs Jul 25, 2024
6dd55a1
update examples
mscroggs Jul 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions examples/lagrange_element.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,9 @@ fn main() {
// Create an array to store the basis function values
let mut basis_values = rlst_dynamic_array4!(f64, element.tabulate_array_shape(0, 1));
// Create array containing the point [1/3, 1/3]
let mut points = rlst_dynamic_array2!(f64, [1, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 1]);
points[[0, 0]] = 1.0 / 3.0;
points[[0, 1]] = 1.0 / 3.0;
points[[1, 0]] = 1.0 / 3.0;
// Tabulate the element's basis functions at the point
element.tabulate(&points, 0, &mut basis_values);
println!(
Expand All @@ -29,7 +29,7 @@ fn main() {

// Set point to [1, 0]
points[[0, 0]] = 1.0;
points[[0, 1]] = 0.0;
points[[1, 0]] = 0.0;
// Tabulate the element's basis functions at the point
element.tabulate(&points, 0, &mut basis_values);
println!(
Expand Down
206 changes: 104 additions & 102 deletions src/ciarlet.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,22 +88,23 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
}

let new_pts = if continuity == Continuity::Discontinuous {
println!("A");
let mut new_pts: EntityPoints<T::Real> = [vec![], vec![], vec![], vec![]];
let mut all_pts = rlst_dynamic_array2![T::Real, [npts, tdim]];
let mut all_pts = rlst_dynamic_array2![T::Real, [tdim, npts]];
for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() {
for _pts in pts_i {
new_pts[i].push(rlst_dynamic_array2![T::Real, [0, tdim]]);
new_pts[i].push(rlst_dynamic_array2![T::Real, [tdim, 0]]);
}
}
let mut row = 0;
let mut col = 0;
for pts_i in interpolation_points.iter() {
for pts in pts_i {
let nrows = pts.shape()[0];
let ncols = pts.shape()[1];
all_pts
.view_mut()
.into_subview([row, 0], [nrows, tdim])
.into_subview([0, col], [tdim, ncols])
.fill_from(pts.view());
row += nrows;
col += ncols;
}
}
new_pts[tdim].push(all_pts);
Expand Down Expand Up @@ -146,8 +147,8 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
let mut dof = 0;
for d in 0..4 {
for (e, pts) in new_pts[d].iter().enumerate() {
if pts.shape()[0] > 0 {
let mut table = rlst_dynamic_array3!(T, [1, pdim, pts.shape()[0]]);
if pts.shape()[1] > 0 {
let mut table = rlst_dynamic_array3!(T, [1, pdim, pts.shape()[1]]);
tabulate_legendre_polynomials(
cell_type,
pts,
Expand Down Expand Up @@ -207,9 +208,9 @@ impl<T: RlstScalar + MatrixInverse> CiarletElement<T> {
let mut dof = 0;
for i in 0..4 {
for pts in &new_pts[i] {
let dofs = (dof..dof + pts.shape()[0]).collect::<Vec<_>>();
let dofs = (dof..dof + pts.shape()[1]).collect::<Vec<_>>();
entity_dofs[i].push(dofs);
dof += pts.shape()[0];
dof += pts.shape()[1];
}
}
let connectivity = reference_cell::connectivity(cell_type);
Expand Down Expand Up @@ -307,7 +308,7 @@ impl<T: RlstScalar + MatrixInverse> FiniteElement for CiarletElement<T> {
);

for d in 0..table.shape()[0] {
for p in 0..points.shape()[0] {
for p in 0..points.shape()[1] {
for j in 0..self.value_size {
for b in 0..self.dim {
// data[d, p, b, j] = inner(self.coefficients[b, j, :], table[d, :, p])
Expand Down Expand Up @@ -384,11 +385,11 @@ mod test {
let e = lagrange::create::<f64>(ReferenceCellType::Interval, 0, Continuity::Discontinuous);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4));
let mut points = rlst_dynamic_array2!(f64, [4, 1]);
let mut points = rlst_dynamic_array2!(f64, [1, 4]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 0.2;
*points.get_mut([2, 0]).unwrap() = 0.4;
*points.get_mut([3, 0]).unwrap() = 1.0;
*points.get_mut([0, 1]).unwrap() = 0.2;
*points.get_mut([0, 2]).unwrap() = 0.4;
*points.get_mut([0, 3]).unwrap() = 1.0;
e.tabulate(&points, 0, &mut data);

for pt in 0..4 {
Expand All @@ -402,21 +403,21 @@ mod test {
let e = lagrange::create::<f64>(ReferenceCellType::Interval, 1, Continuity::Standard);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 4));
let mut points = rlst_dynamic_array2!(f64, [4, 1]);
let mut points = rlst_dynamic_array2!(f64, [1, 4]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 0.2;
*points.get_mut([2, 0]).unwrap() = 0.4;
*points.get_mut([3, 0]).unwrap() = 1.0;
*points.get_mut([0, 1]).unwrap() = 0.2;
*points.get_mut([0, 2]).unwrap() = 0.4;
*points.get_mut([0, 3]).unwrap() = 1.0;
e.tabulate(&points, 0, &mut data);

for pt in 0..4 {
assert_relative_eq!(
*data.get([0, pt, 0, 0]).unwrap(),
1.0 - *points.get([pt, 0]).unwrap()
1.0 - *points.get([0, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 1, 0]).unwrap(),
*points.get([pt, 0]).unwrap()
*points.get([0, pt]).unwrap()
);
}
check_dofs(e);
Expand All @@ -428,19 +429,19 @@ mod test {
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));

let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 0.5;
*points.get_mut([3, 1]).unwrap() = 0.0;
*points.get_mut([4, 0]).unwrap() = 0.0;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.5;
*points.get_mut([5, 1]).unwrap() = 0.5;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 0.5;
*points.get_mut([1, 3]).unwrap() = 0.0;
*points.get_mut([0, 4]).unwrap() = 0.0;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.5;
*points.get_mut([1, 5]).unwrap() = 0.5;

e.tabulate(&points, 0, &mut data);

Expand All @@ -455,33 +456,33 @@ mod test {
let e = lagrange::create::<f64>(ReferenceCellType::Triangle, 1, Continuity::Standard);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 0.5;
*points.get_mut([3, 1]).unwrap() = 0.0;
*points.get_mut([4, 0]).unwrap() = 0.0;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.5;
*points.get_mut([5, 1]).unwrap() = 0.5;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 0.5;
*points.get_mut([1, 3]).unwrap() = 0.0;
*points.get_mut([0, 4]).unwrap() = 0.0;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.5;
*points.get_mut([1, 5]).unwrap() = 0.5;
e.tabulate(&points, 0, &mut data);

for pt in 0..6 {
assert_relative_eq!(
*data.get([0, pt, 0, 0]).unwrap(),
1.0 - *points.get([pt, 0]).unwrap() - *points.get([pt, 1]).unwrap()
1.0 - *points.get([0, pt]).unwrap() - *points.get([1, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 1, 0]).unwrap(),
*points.get([pt, 0]).unwrap()
*points.get([0, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 2, 0]).unwrap(),
*points.get([pt, 1]).unwrap()
*points.get([1, pt]).unwrap()
);
}
check_dofs(e);
Expand Down Expand Up @@ -551,19 +552,19 @@ mod test {
);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 0.5;
*points.get_mut([3, 1]).unwrap() = 0.0;
*points.get_mut([4, 0]).unwrap() = 0.0;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.5;
*points.get_mut([5, 1]).unwrap() = 0.5;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 0.5;
*points.get_mut([1, 3]).unwrap() = 0.0;
*points.get_mut([0, 4]).unwrap() = 0.0;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.5;
*points.get_mut([1, 5]).unwrap() = 0.5;
e.tabulate(&points, 0, &mut data);

for pt in 0..6 {
Expand All @@ -577,38 +578,38 @@ mod test {
let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 1, Continuity::Standard);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 1.0;
*points.get_mut([3, 1]).unwrap() = 1.0;
*points.get_mut([4, 0]).unwrap() = 0.25;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.3;
*points.get_mut([5, 1]).unwrap() = 0.2;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 1.0;
*points.get_mut([1, 3]).unwrap() = 1.0;
*points.get_mut([0, 4]).unwrap() = 0.25;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.3;
*points.get_mut([1, 5]).unwrap() = 0.2;

e.tabulate(&points, 0, &mut data);

for pt in 0..6 {
assert_relative_eq!(
*data.get([0, pt, 0, 0]).unwrap(),
(1.0 - *points.get([pt, 0]).unwrap()) * (1.0 - *points.get([pt, 1]).unwrap())
(1.0 - *points.get([0, pt]).unwrap()) * (1.0 - *points.get([1, pt]).unwrap())
);
assert_relative_eq!(
*data.get([0, pt, 1, 0]).unwrap(),
*points.get([pt, 0]).unwrap() * (1.0 - *points.get([pt, 1]).unwrap())
*points.get([0, pt]).unwrap() * (1.0 - *points.get([1, pt]).unwrap())
);
assert_relative_eq!(
*data.get([0, pt, 2, 0]).unwrap(),
(1.0 - *points.get([pt, 0]).unwrap()) * *points.get([pt, 1]).unwrap()
(1.0 - *points.get([0, pt]).unwrap()) * *points.get([1, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 3, 0]).unwrap(),
*points.get([pt, 0]).unwrap() * *points.get([pt, 1]).unwrap()
*points.get([0, pt]).unwrap() * *points.get([1, pt]).unwrap()
);
}
check_dofs(e);
Expand All @@ -619,24 +620,25 @@ mod test {
let e = lagrange::create::<f64>(ReferenceCellType::Quadrilateral, 2, Continuity::Standard);
assert_eq!(e.value_size(), 1);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 1.0;
*points.get_mut([3, 1]).unwrap() = 1.0;
*points.get_mut([4, 0]).unwrap() = 0.25;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.3;
*points.get_mut([5, 1]).unwrap() = 0.2;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 1.0;
*points.get_mut([1, 3]).unwrap() = 1.0;
*points.get_mut([0, 4]).unwrap() = 0.25;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.3;
*points.get_mut([1, 5]).unwrap() = 0.2;

e.tabulate(&points, 0, &mut data);

for pt in 0..6 {
let x = *points.get([pt, 0]).unwrap();
let y = *points.get([pt, 1]).unwrap();
let x = *points.get([0, pt]).unwrap();
let y = *points.get([1, pt]).unwrap();
assert_relative_eq!(
*data.get([0, pt, 0, 0]).unwrap(),
(1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y),
Expand Down Expand Up @@ -691,45 +693,45 @@ mod test {
let e = raviart_thomas::create(ReferenceCellType::Triangle, 1, Continuity::Standard);
assert_eq!(e.value_size(), 2);
let mut data = rlst_dynamic_array4!(f64, e.tabulate_array_shape(0, 6));
let mut points = rlst_dynamic_array2!(f64, [6, 2]);
let mut points = rlst_dynamic_array2!(f64, [2, 6]);
*points.get_mut([0, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 0.0;
*points.get_mut([1, 0]).unwrap() = 1.0;
*points.get_mut([1, 0]).unwrap() = 0.0;
*points.get_mut([0, 1]).unwrap() = 1.0;
*points.get_mut([1, 1]).unwrap() = 0.0;
*points.get_mut([2, 0]).unwrap() = 0.0;
*points.get_mut([2, 1]).unwrap() = 1.0;
*points.get_mut([3, 0]).unwrap() = 0.5;
*points.get_mut([3, 1]).unwrap() = 0.0;
*points.get_mut([4, 0]).unwrap() = 0.0;
*points.get_mut([4, 1]).unwrap() = 0.5;
*points.get_mut([5, 0]).unwrap() = 0.5;
*points.get_mut([5, 1]).unwrap() = 0.5;
*points.get_mut([0, 2]).unwrap() = 0.0;
*points.get_mut([1, 2]).unwrap() = 1.0;
*points.get_mut([0, 3]).unwrap() = 0.5;
*points.get_mut([1, 3]).unwrap() = 0.0;
*points.get_mut([0, 4]).unwrap() = 0.0;
*points.get_mut([1, 4]).unwrap() = 0.5;
*points.get_mut([0, 5]).unwrap() = 0.5;
*points.get_mut([1, 5]).unwrap() = 0.5;
e.tabulate(&points, 0, &mut data);

for pt in 0..6 {
assert_relative_eq!(
*data.get([0, pt, 0, 0]).unwrap(),
-*points.get([pt, 0]).unwrap()
-*points.get([0, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 0, 1]).unwrap(),
-*points.get([pt, 1]).unwrap()
-*points.get([1, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 1, 0]).unwrap(),
*points.get([pt, 0]).unwrap() - 1.0
*points.get([0, pt]).unwrap() - 1.0
);
assert_relative_eq!(
*data.get([0, pt, 1, 1]).unwrap(),
*points.get([pt, 1]).unwrap()
*points.get([1, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 2, 0]).unwrap(),
-*points.get([pt, 0]).unwrap()
-*points.get([0, pt]).unwrap()
);
assert_relative_eq!(
*data.get([0, pt, 2, 1]).unwrap(),
1.0 - *points.get([pt, 1]).unwrap()
1.0 - *points.get([1, pt]).unwrap()
);
}
check_dofs(e);
Expand Down
Loading
Loading