Skip to content

Commit

Permalink
get element working without Array2D
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Oct 17, 2023
1 parent fb122ce commit 4b5affa
Show file tree
Hide file tree
Showing 12 changed files with 166 additions and 361 deletions.
23 changes: 16 additions & 7 deletions bem/src/assembly/batched.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,16 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz
let gdim = test_grid.geometry().dim();

// TODO: move this to grid.get_compute_points_function(test_points)
let test_compute_points = |cell: usize, pts: &mut rlst_dense::Matrix<f64, rlst_dense::base_matrix::BaseMatrix<f64, rlst_dense::VectorContainer<f64>, rlst_dense::Dynamic>, rlst_dense::Dynamic>| {
let test_compute_points = |cell: usize,
pts: &mut rlst_dense::Matrix<
f64,
rlst_dense::base_matrix::BaseMatrix<
f64,
rlst_dense::VectorContainer<f64>,
rlst_dense::Dynamic,
>,
rlst_dense::Dynamic,
>| {
for p in 0..NPTS_TEST {
for i in 0..gdim {
unsafe {
Expand Down Expand Up @@ -454,8 +463,8 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize>
test_start = test_end
}

let numthreads = test_cells.len();
let r: usize = (0..numthreads)
let numtasks = test_cells.len();
let r: usize = (0..numtasks)
.into_par_iter()
.map(&|t| {
assemble_batch_nonadjacent::<NPTS_TEST, NPTS_TRIAL>(
Expand All @@ -476,7 +485,7 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize>
)
})
.sum();
assert_eq!(r, numthreads);
assert_eq!(r, numtasks);
}
}
}
Expand Down Expand Up @@ -616,8 +625,8 @@ pub fn assemble_singular<'a>(
start = end;
}

let numthreads = cell_blocks.len();
let r: usize = (0..numthreads)
let numtasks = cell_blocks.len();
let r: usize = (0..numtasks)
.into_par_iter()
.map(&|t| {
assemble_batch_singular(
Expand All @@ -636,7 +645,7 @@ pub fn assemble_singular<'a>(
)
})
.sum();
assert_eq!(r, numthreads);
assert_eq!(r, numtasks);
}
}
}
Expand Down
18 changes: 9 additions & 9 deletions element/src/element.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,13 @@

use crate::cell::create_cell;
use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials};
use bempp_tools::arrays::{AdjacencyList, Array2D, Array3D};
use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess};
use bempp_tools::arrays::{AdjacencyList, Array3D, Mat};
use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess};
use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{Continuity, ElementFamily, FiniteElement, MapType};
use rlst_algorithms::linalg::LinAlg;
use rlst_algorithms::traits::inverse::Inverse;
use rlst_dense::{RandomAccessByRef, RandomAccessMut};
use rlst_dense::{rlst_dynamic_mat, RandomAccessByRef, RandomAccessMut, Shape};
pub mod lagrange;
pub mod raviart_thomas;

Expand Down Expand Up @@ -37,7 +37,7 @@ impl CiarletElement {
degree: usize,
value_shape: Vec<usize>,
polynomial_coeffs: Array3D<f64>,
interpolation_points: [Vec<Array2D<f64>>; 4],
interpolation_points: [Vec<Mat<f64>>; 4],
interpolation_weights: [Vec<Array3D<f64>>; 4],
map_type: MapType,
continuity: Continuity,
Expand Down Expand Up @@ -69,12 +69,12 @@ impl CiarletElement {
}

let new_pts = if continuity == Continuity::Discontinuous {
let mut new_pts = [vec![], vec![], vec![], vec![]];
let mut new_pts: [Vec<Mat<f64>>; 4] = [vec![], vec![], vec![], vec![]];
let mut pn = 0;
let mut all_pts = Array2D::<f64>::new((npts, tdim));
let mut all_pts = rlst_dynamic_mat![f64, (npts, tdim)];
for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() {
for _pts in pts_i {
new_pts[i].push(Array2D::<f64>::new((0, tdim)));
new_pts[i].push(rlst_dynamic_mat![f64, (0, tdim)]);
}
}
for pts_i in interpolation_points.iter() {
Expand Down Expand Up @@ -240,9 +240,9 @@ impl FiniteElement for CiarletElement {
fn dim(&self) -> usize {
self.dim
}
fn tabulate<'a>(
fn tabulate<T: RandomAccessByRef<Item = f64> + Shape>(
&self,
points: &impl Array2DAccess<'a, f64>,
points: &T,
nderivs: usize,
data: &mut impl Array4DAccess<f64>,
) {
Expand Down
41 changes: 21 additions & 20 deletions element/src/element/lagrange.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

use crate::element::{create_cell, CiarletElement};
use crate::polynomials::polynomial_count;
use bempp_tools::arrays::{Array2D, Array3D};
use bempp_tools::arrays::{to_matrix, Array3D};
use bempp_traits::arrays::Array3DAccess;
use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{Continuity, ElementFamily, MapType};
Expand All @@ -29,18 +29,18 @@ pub fn create(
}
for d in 0..tdim {
for _e in 0..cell.entity_count(d) {
x[d].push(Array2D::<f64>::new((0, tdim)));
x[d].push(to_matrix(&[], (0, tdim)));
m[d].push(Array3D::<f64>::new((0, 1, 0)));
}
}
x[tdim].push(Array2D::<f64>::from_data(cell.midpoint(), (1, tdim)));
x[tdim].push(to_matrix(&cell.midpoint(), (1, tdim)));
m[tdim].push(Array3D::<f64>::from_data(vec![1.0], (1, 1, 1)));
} else {
// TODO: GLL points
for e in 0..cell.entity_count(0) {
let mut pts = vec![0.0; tdim];
pts.copy_from_slice(&cell.vertices()[e * tdim..(e + 1) * tdim]);
x[0].push(Array2D::<f64>::from_data(pts, (1, tdim)));
x[0].push(to_matrix(&pts, (1, tdim)));
m[0].push(Array3D::<f64>::from_data(vec![1.0], (1, 1, 1)));
}
for e in 0..cell.entity_count(1) {
Expand All @@ -56,7 +56,7 @@ pub fn create(
pts[(i - 1) * tdim + j] = v0[j] + i as f64 / degree as f64 * (v1[j] - v0[j]);
}
}
x[1].push(Array2D::<f64>::from_data(pts, (degree - 1, tdim)));
x[1].push(to_matrix(&pts, (degree - 1, tdim)));
m[1].push(Array3D::<f64>::from_data(
ident,
(degree - 1, 1, degree - 1),
Expand Down Expand Up @@ -119,7 +119,7 @@ pub fn create(
for i in 0..npts {
ident[i * npts + i] = 1.0;
}
x[2].push(Array2D::<f64>::from_data(pts, (npts, tdim)));
x[2].push(to_matrix(&pts, (npts, tdim)));
m[2].push(Array3D::<f64>::from_data(ident, (npts, 1, npts)));
start += nvertices;
}
Expand All @@ -143,9 +143,10 @@ mod test {
use crate::cell::*;
use crate::element::lagrange::*;
use approx::*;
use bempp_tools::arrays::{Array2D, Array4D};
use bempp_traits::arrays::{Array2DAccess, Array4DAccess};
use bempp_tools::arrays::Array4D;
use bempp_traits::arrays::Array4DAccess;
use bempp_traits::element::FiniteElement;
use rlst_dense::RandomAccessByRef;

fn check_dofs(e: impl FiniteElement) {
let cell_dim = match e.cell_type() {
Expand Down Expand Up @@ -182,7 +183,7 @@ mod test {
let e = create(ReferenceCellType::Interval, 0, Continuity::Discontinuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 4));
let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1));
let points = to_matrix(&vec![0.0, 0.2, 0.4, 1.0], (4, 1));
e.tabulate(&points, 0, &mut data);

for pt in 0..4 {
Expand All @@ -196,7 +197,7 @@ mod test {
let e = create(ReferenceCellType::Interval, 1, Continuity::Continuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 4));
let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1));
let points = to_matrix(&vec![0.0, 0.2, 0.4, 1.0], (4, 1));
e.tabulate(&points, 0, &mut data);

for pt in 0..4 {
Expand All @@ -214,8 +215,8 @@ mod test {
let e = create(ReferenceCellType::Triangle, 0, Continuity::Discontinuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand All @@ -231,8 +232,8 @@ mod test {
let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand Down Expand Up @@ -312,8 +313,8 @@ mod test {
);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand All @@ -329,8 +330,8 @@ mod test {
let e = create(ReferenceCellType::Quadrilateral, 1, Continuity::Continuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand Down Expand Up @@ -361,8 +362,8 @@ mod test {
let e = create(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous);
assert_eq!(e.value_size(), 1);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand Down
17 changes: 9 additions & 8 deletions element/src/element/raviart_thomas.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

use crate::element::{create_cell, CiarletElement};
use crate::polynomials::polynomial_count;
use bempp_tools::arrays::{Array2D, Array3D};
use bempp_tools::arrays::{to_matrix, Array3D};
use bempp_traits::arrays::Array3DAccess;
use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{Continuity, ElementFamily, MapType};
Expand Down Expand Up @@ -45,7 +45,7 @@ pub fn create(
let mut x = [vec![], vec![], vec![], vec![]];
let mut m = [vec![], vec![], vec![], vec![]];
for _e in 0..cell.entity_count(0) {
x[0].push(Array2D::<f64>::new((0, tdim)));
x[0].push(to_matrix(&[], (0, tdim)));
m[0].push(Array3D::<f64>::new((0, 2, 0)));
}

Expand All @@ -61,12 +61,12 @@ pub fn create(
}
mat[0] = v0[1] - v1[1];
mat[1] = v1[0] - v0[0];
x[1].push(Array2D::<f64>::from_data(pts, (1, tdim)));
x[1].push(to_matrix(&pts, (1, tdim)));
m[1].push(Array3D::<f64>::from_data(mat, (1, 2, 1)));
}

for _e in 0..cell.entity_count(2) {
x[2].push(Array2D::<f64>::new((0, tdim)));
x[2].push(to_matrix(&[], (0, tdim)));
m[2].push(Array3D::<f64>::new((0, 2, 0)));
}

Expand All @@ -89,9 +89,10 @@ mod test {
use crate::cell::*;
use crate::element::raviart_thomas::*;
use approx::*;
use bempp_tools::arrays::{Array2D, Array4D};
use bempp_traits::arrays::{Array2DAccess, Array4DAccess};
use bempp_tools::arrays::Array4D;
use bempp_traits::arrays::Array4DAccess;
use bempp_traits::element::FiniteElement;
use rlst_dense::RandomAccessByRef;

fn check_dofs(e: impl FiniteElement) {
let cell_dim = match e.cell_type() {
Expand Down Expand Up @@ -128,8 +129,8 @@ mod test {
let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous);
assert_eq!(e.value_size(), 2);
let mut data = Array4D::<f64>::new(e.tabulate_array_shape(0, 6));
let points = Array2D::from_data(
vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
let points = to_matrix(
&vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5],
(6, 2),
);
e.tabulate(&points, 0, &mut data);
Expand Down
Loading

0 comments on commit 4b5affa

Please sign in to comment.