diff --git a/Cargo.toml b/Cargo.toml index 296dab4..6dbc499 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,3 +15,4 @@ rlst = { git = "https://github.com/linalg-rs/rlst.git" } rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } rlst-proc-macro = { git = "https://github.com/linalg-rs/rlst.git" } rlst-common = { git = "https://github.com/linalg-rs/rlst.git" } +log = "0.4" diff --git a/src/grid/surface_triangle_grid.rs b/src/grid/surface_triangle_grid.rs index 5688791..e584adf 100644 --- a/src/grid/surface_triangle_grid.rs +++ b/src/grid/surface_triangle_grid.rs @@ -6,258 +6,3 @@ pub mod grid; pub mod reference_map; pub mod topology; pub mod vertex; - -// use std::{collections::HashMap, iter::Copied, marker::PhantomData, ops::Add}; - -// use num::{traits::Float, One}; -// use rlst::{rlst_dynamic_array2, rlst_static_array, Array}; -// use rlst_common::types::Scalar; -// use rlst_dense::{ -// array::DynamicArray, -// rlst_array_from_slice1, rlst_array_from_slice_mut1, rlst_array_from_slice_mut2, -// traits::{DefaultIteratorMut, RandomAccessMut, RawAccess}, -// }; -// use rlst_proc_macro::rlst_static_type; - -// use crate::traits::{CellType, GridType, TopologyType, VertexType}; - -// pub struct TriangleVertex<'a, T: Float + Scalar> { -// id: usize, -// index: usize, -// data: &'a [T; 3], -// } - -// impl<'a, T: Float + Scalar> VertexType for TriangleVertex<'a, T> { -// type T = T; - -// fn coords(&self, data: &mut [Self::T]) { -// for (out_data, &in_data) in data.iter_mut().zip(self.data) { -// *out_data = in_data; -// } -// } - -// fn index(&self) -> usize { -// self.index -// } - -// fn id(&self) -> usize { -// self.id -// } -// } - -// pub struct TriangleCell<'a, T: Float + Scalar> { -// id: usize, -// index: usize, -// data: &'a [usize; 3], -// _marker: PhantomData, -// } - -// impl<'a, T: Float + Scalar> CellType for TriangleCell<'a, T> { -// type Topology<'top> = TriangleTopology<'top, T> where Self: 'top; - -// fn id(&self) -> usize { -// self.id -// } - -// fn index(&self) -> usize { -// self.index -// } - -// fn topology(&self) -> Self::Topology<'_> { -// TriangleTopology::new(self) -// } -// } - -// pub struct TriangleTopology<'a, T: Float + Scalar> { -// cell: &'a TriangleCell<'a, T>, -// _marker: PhantomData, -// } - -// impl<'a, T: Float + Scalar> TriangleTopology<'a, T> { -// pub fn new(cell: &'a TriangleCell<'a, T>) -> Self { -// Self { -// cell, -// _marker: PhantomData, -// } -// } -// } - -// impl<'a, T: Float + Scalar> TopologyType for TriangleTopology<'a, T> { -// type Grid = TriangleSurfaceGrid; - -// type VertexIndexIter<'iter> = Copied> where Self: 'iter; - -// fn vertex_indices(&self) -> Self::VertexIndexIter<'_> { -// self.cell.data.iter().copied() -// } -// } - -// // pub struct TriangleMap { -// // v0: rlst_static_type!(T, 3), -// // gradient: rlst_static_type!(T, 3, 2), -// // } - -// // impl TriangleMap { -// // pub fn new(cell: &TriangleCell, grid: &TriangleSurfaceGrid) -> Self { -// // let mut v0 = rlst_static_array!(T, 3); -// // v0.fill_from(grid.vertices[cell.vertices()[0]].data().view()); - -// // let mut gradient = rlst_static_array!(T, 3, 2); -// // let vertices = [ -// // &grid.vertices[cell.vertices()[1]], -// // &grid.vertices[cell.vertices()[2]], -// // ]; - -// // for (mut gradient_col, vertex) in gradient.col_iter_mut().zip(vertices) { -// // gradient_col.fill_from( -// // vertex.data().view().add( -// // v0.view() -// // .add(vertex.data().view().scalar_mul(-::one())), -// // ), -// // ); -// // } - -// // Self { v0, gradient } -// // } -// // } - -// pub struct TriangleSurfaceGrid { -// vertices: Vec<[T; 3]>, -// cells: Vec<[usize; 3]>, -// ids_from_vertex_indices: Vec, -// ids_from_cell_indices: Vec, -// vertex_ids: HashMap, -// cell_ids: HashMap, -// } - -// impl Default for TriangleSurfaceGrid { -// fn default() -> Self { -// Self::new() -// } -// } - -// impl TriangleSurfaceGrid { -// pub fn new() -> Self { -// Self { -// vertices: Vec::new(), -// cells: Vec::new(), -// ids_from_vertex_indices: Vec::new(), -// ids_from_cell_indices: Vec::new(), -// vertex_ids: HashMap::new(), -// cell_ids: HashMap::new(), -// } -// } - -// pub fn new_with_capacity(nvertices: usize, ncells: usize) -> Self { -// Self { -// vertices: Vec::with_capacity(nvertices), -// cells: Vec::with_capacity(ncells), -// ids_from_vertex_indices: Vec::with_capacity(nvertices), -// ids_from_cell_indices: Vec::with_capacity(ncells), -// vertex_ids: HashMap::new(), -// cell_ids: HashMap::new(), -// } -// } - -// pub fn add_vertex(&mut self, id: usize, data: [T; 3]) { -// let index = self.vertices.len(); -// self.vertices.push(data); -// self.ids_from_vertex_indices.push(id); -// self.vertex_ids.insert(id, index); -// } - -// pub fn add_cell(&mut self, id: usize, vertex_ids: [usize; 3]) { -// let index = self.cells.len(); -// self.cells.push([ -// self.vertex_ids[&vertex_ids[0]], -// self.vertex_ids[&vertex_ids[1]], -// self.vertex_ids[&vertex_ids[2]], -// ]); -// self.ids_from_cell_indices.push(id); -// self.cell_ids.insert(id, index); -// } -// } - -// impl GridType for TriangleSurfaceGrid { -// type T = T; - -// type Vertex<'a> = TriangleVertex<'a, T> where Self: 'a; -// type Cell<'a> = TriangleCell<'a, T> where Self: 'a; - -// type Edge = [usize; 2]; - -// type Face = (); - -// fn number_of_vertices(&self) -> usize { -// self.vertices.len() -// } -// fn number_of_cells(&self) -> usize { -// self.cells.len() -// } - -// fn vertex_index_from_id(&self, id: usize) -> usize { -// self.vertex_ids[&id] -// } -// fn vertex_id_from_index(&self, index: usize) -> usize { -// self.ids_from_vertex_indices[index] -// } - -// fn cell_index_from_id(&self, id: usize) -> usize { -// self.cell_ids[&id] -// } -// fn cell_id_from_index(&self, index: usize) -> usize { -// self.ids_from_cell_indices[index] -// } - -// fn vertex_from_index(&self, index: usize) -> Self::Vertex<'_> { -// TriangleVertex { -// id: self.vertex_id_from_index(index), -// index, -// data: &self.vertices[index], -// } -// } - -// fn cell_from_index(&self, index: usize) -> Self::Cell<'_> { -// TriangleCell { -// id: self.cell_id_from_index(index), -// index, -// data: &self.cells[index], -// _marker: PhantomData, -// } -// } -// } - -// #[cfg(test)] -// mod test { -// use crate::traits::{CellType, GridType, TopologyType}; - -// use super::TriangleSurfaceGrid; - -// #[test] -// fn test_grid() { -// let mut grid = TriangleSurfaceGrid::::new(); - -// grid.add_vertex(0, [0.0, 0.0, 0.0]); -// grid.add_vertex(1, [1.0, 0.0, 0.0]); -// grid.add_vertex(2, [0.0, 1.0, 0.0]); -// grid.add_vertex(3, [1.0, 1.0, 0.0]); - -// grid.add_cell(0, [0, 1, 2]); -// grid.add_cell(0, [1, 3, 2]); - -// for vertex in grid.iter_all_vertices() { -// println!("{:#?}", vertex.data) -// } - -// for cell in grid.iter_all_cells() { -// for (local_index, vertex_index) in cell.topology().vertex_indices().enumerate() { -// println!( -// "Cell: {}, Vertex: {}, {}", -// cell.index(), -// local_index, -// vertex_index -// ) -// } -// } -// } -// } diff --git a/src/grid/surface_triangle_grid/geometry.rs b/src/grid/surface_triangle_grid/geometry.rs index 1de76ee..bcc0983 100644 --- a/src/grid/surface_triangle_grid/geometry.rs +++ b/src/grid/surface_triangle_grid/geometry.rs @@ -8,7 +8,7 @@ use rlst_dense::traits::DefaultIterator; use crate::{ traits::{GeometryType, GridType}, - types::vertex_iterator::PointIterator, + types::point_iterator::PointIterator, }; use super::grid::TriangleSurfaceGrid; diff --git a/src/grid/surface_triangle_grid/grid.rs b/src/grid/surface_triangle_grid/grid.rs index e56e6a5..f1a3a1d 100644 --- a/src/grid/surface_triangle_grid/grid.rs +++ b/src/grid/surface_triangle_grid/grid.rs @@ -10,11 +10,7 @@ use rlst_proc_macro::rlst_static_type; use crate::traits::*; -use super::{ - cell::TriangleCell, - reference_map::{TriangleReferenceMap, TriangleReferenceMapIterator}, - vertex::TriangleVertex, -}; +use super::{cell::TriangleCell, reference_map::TriangleReferenceMap, vertex::TriangleVertex}; pub struct TriangleSurfaceGrid { pub(crate) vertices: Vec<[T; 3]>, @@ -23,8 +19,8 @@ pub struct TriangleSurfaceGrid { ids_from_cell_indices: Vec, vertex_ids: HashMap, cell_ids: HashMap, - edge_to_cells: Vec>, - point_to_cells: Vec>, + edge_to_cells: Vec>>, + point_to_cells: Vec>>, pub(crate) cell_to_edges: Vec<[usize; 3]>, pub(crate) jacobians: Vec, pub(crate) volumes: Vec, @@ -108,7 +104,7 @@ impl TriangleSurfaceGrid { .resize_with(self.cells.len(), Default::default); let edge_local: [(usize, usize); 3] = [(1, 2), (0, 2), (0, 1)]; let mut edge_connectivity = - HashMap::<(usize, usize), (usize, Vec)>::new(); + HashMap::<(usize, usize), (usize, Vec>)>::new(); for (cell_index, cell_vertices) in self.cells.iter().enumerate() { for (local_index, &(first, second)) in edge_local.iter().enumerate() { let mut first = cell_vertices[first]; @@ -119,11 +115,11 @@ impl TriangleSurfaceGrid { if let Some((edge_index, adjacent_cells)) = edge_connectivity.get_mut(&(first, second)) { - adjacent_cells.push(CellLocalIndexPair::new(cell_index, local_index)); + adjacent_cells.push(CellLocalIndexPair::::new(cell_index, local_index)); self.cell_to_edges[cell_index][local_index] = *edge_index; } else { - let mut adjacent_cells = Vec::::with_capacity(2); - adjacent_cells.push(CellLocalIndexPair::new(cell_index, local_index)); + let mut adjacent_cells = Vec::>::with_capacity(2); + adjacent_cells.push(CellLocalIndexPair::::new(cell_index, local_index)); self.cell_to_edges[cell_index][local_index] = nedges; edge_connectivity.insert((first, second), (nedges, adjacent_cells)); nedges += 1; @@ -143,7 +139,7 @@ impl TriangleSurfaceGrid { for (cell_index, cell) in self.cells.iter().enumerate() { for (local_index, &vertex_index) in cell.iter().enumerate() { self.point_to_cells[vertex_index] - .push(CellLocalIndexPair::new(cell_index, local_index)); + .push(CellLocalIndexPair::::new(cell_index, local_index)); } } } @@ -200,7 +196,7 @@ impl TriangleSurfaceGrid { impl GridType for TriangleSurfaceGrid { type T = T; - + type IndexType = usize; type Point<'a> = TriangleVertex<'a, T> where Self: 'a; type Cell<'a> = TriangleCell<'a, T> where Self: 'a; @@ -208,11 +204,6 @@ impl GridType for TriangleSurfaceGrid { where Self: 'a; - type ReferenceMapIterator<'a, Iter: std::iter::Iterator> = TriangleReferenceMapIterator<'a, Self, Iter> - where - Self: 'a, - Iter: 'a; - fn number_of_vertices(&self) -> usize { self.vertices.len() } @@ -254,31 +245,19 @@ impl GridType for TriangleSurfaceGrid { fn reference_to_physical_map<'a>( &'a self, reference_points: &'a [Self::T], - cell_index: usize, ) -> Self::ReferenceMap<'a> { - TriangleReferenceMap::new(reference_points, cell_index, self) - } - - fn iter_reference_to_physical_map<'a, Iter: std::iter::Iterator + 'a>( - &'a self, - reference_points: &'a [Self::T], - iter: Iter, - ) -> Self::ReferenceMapIterator<'a, Iter> - where - Self: 'a, - { - TriangleReferenceMapIterator::new(iter, reference_points, self) + TriangleReferenceMap::new(reference_points, self) } - fn edge_to_cells(&self, edge_index: usize) -> &[CellLocalIndexPair] { + fn edge_to_cells(&self, edge_index: usize) -> &[CellLocalIndexPair] { self.edge_to_cells[edge_index].as_slice() } - fn point_to_cells(&self, point_index: usize) -> &[CellLocalIndexPair] { - self.point_to_cells[point_index].as_slice() + fn vertex_to_cells(&self, vertex_index: usize) -> &[CellLocalIndexPair] { + self.point_to_cells[vertex_index].as_slice() } - fn face_to_cells(&self, _face_index: usize) -> &[CellLocalIndexPair] { + fn face_to_cells(&self, _face_index: usize) -> &[CellLocalIndexPair] { std::unimplemented!() } } @@ -341,12 +320,12 @@ mod test { reference_points[[0, 1]] = 1.0; reference_points[[1, 1]] = 0.0; - for map in - grid.iter_reference_to_physical_map(reference_points.data(), 0..grid.number_of_cells()) - { + let map = grid.reference_to_physical_map(reference_points.data()); + for cell_index in 0..grid.number_of_cells() { let mut points = rlst_dynamic_array2!(f64, [3, map.number_of_reference_points()]); for point_index in 0..map.number_of_reference_points() { map.reference_to_physical( + cell_index, point_index, points.view_mut().slice(1, point_index).data_mut(), ); diff --git a/src/grid/surface_triangle_grid/reference_map.rs b/src/grid/surface_triangle_grid/reference_map.rs index 0d2ed83..7c759d1 100644 --- a/src/grid/surface_triangle_grid/reference_map.rs +++ b/src/grid/surface_triangle_grid/reference_map.rs @@ -14,21 +14,15 @@ use super::grid::TriangleSurfaceGrid; pub struct TriangleReferenceMap<'a, T: Float + Scalar> { reference_points: SliceArray<'a, T, 2>, - cell_index: usize, grid: &'a TriangleSurfaceGrid, } impl<'a, T: Float + Scalar> TriangleReferenceMap<'a, T> { - pub fn new( - reference_points: &'a [T], - cell_index: usize, - grid: &'a TriangleSurfaceGrid, - ) -> Self { + pub fn new(reference_points: &'a [T], grid: &'a TriangleSurfaceGrid) -> Self { assert_eq!(reference_points.len() % 2, 0); let npoints = reference_points.len() / 2; Self { reference_points: rlst_array_from_slice2!(T, reference_points, [2, npoints]), - cell_index, grid, } } @@ -49,10 +43,15 @@ impl<'a, T: Float + Scalar> ReferenceMapType for TriangleReferenceMap<'a, T> { self.reference_points.shape()[1] } - fn reference_to_physical(&self, point_index: usize, value: &mut [::T]) { + fn reference_to_physical( + &self, + cell_index: usize, + point_index: usize, + value: &mut [::T], + ) { assert_eq!(value.len(), 3); - let v0 = self.grid.vertices[self.grid.cells[self.cell_index][0]]; - let jacobian = &self.grid.jacobians[self.cell_index]; + let v0 = self.grid.vertices[self.grid.cells[cell_index][0]]; + let jacobian = &self.grid.jacobians[cell_index]; for (index, val_out) in value.iter_mut().enumerate() { *val_out = v0[index] + jacobian[[index, 0]] * self.reference_points[[0, point_index]] @@ -60,57 +59,25 @@ impl<'a, T: Float + Scalar> ReferenceMapType for TriangleReferenceMap<'a, T> { } } - fn jacobian(&self, _point_index: usize, value: &mut [::T]) { - for (elem_in, val_out) in self.grid.jacobians[self.cell_index] - .iter() - .zip(value.iter_mut()) - { + fn jacobian( + &self, + cell_index: usize, + _point_index: usize, + value: &mut [::T], + ) { + for (elem_in, val_out) in self.grid.jacobians[cell_index].iter().zip(value.iter_mut()) { *val_out = elem_in } } - fn normal(&self, _point_index: usize, value: &mut [::T]) { - for (elem_in, val_out) in self.grid.normals[self.cell_index] - .iter() - .zip(value.iter_mut()) - { + fn normal( + &self, + cell_index: usize, + _point_index: usize, + value: &mut [::T], + ) { + for (elem_in, val_out) in self.grid.normals[cell_index].iter().zip(value.iter_mut()) { *val_out = elem_in } } } - -pub struct TriangleReferenceMapIterator<'a, Grid: GridType, Iter: std::iter::Iterator> -{ - iter: Iter, - reference_points: &'a [Grid::T], - grid: &'a Grid, -} - -impl<'a, Grid: GridType, Iter: std::iter::Iterator> - TriangleReferenceMapIterator<'a, Grid, Iter> -{ - pub fn new(iter: Iter, reference_points: &'a [Grid::T], grid: &'a Grid) -> Self { - Self { - iter, - reference_points, - grid, - } - } -} - -impl<'a, Grid: GridType, Iter: std::iter::Iterator> Iterator - for TriangleReferenceMapIterator<'a, Grid, Iter> -{ - type Item = Grid::ReferenceMap<'a>; - - fn next(&mut self) -> Option { - if let Some(cell_index) = self.iter.next() { - Some( - self.grid - .reference_to_physical_map(self.reference_points, cell_index), - ) - } else { - None - } - } -} diff --git a/src/grid_impl.rs b/src/grid_impl.rs index 0816410..bad5ce7 100644 --- a/src/grid_impl.rs +++ b/src/grid_impl.rs @@ -1,5 +1,7 @@ //! Grid implementation +pub mod common; +pub mod flat_triangle_grid; pub mod grid; pub mod mixed_grid; pub mod single_element_grid; diff --git a/src/grid_impl/common.rs b/src/grid_impl/common.rs new file mode 100644 index 0000000..6d868e4 --- /dev/null +++ b/src/grid_impl/common.rs @@ -0,0 +1,66 @@ +//! Functionality common to multiple grid implementations + +use crate::grid_impl::traits::Geometry; +use rlst_common::types::Scalar; +use rlst_dense::traits::{Shape, UnsafeRandomAccessByRef}; + +pub fn compute_point + Shape<4>>( + geometry: &impl Geometry, + table: Table, + cell_index: usize, + point_index: usize, + point: &mut [T], +) { + assert_eq!(geometry.dim(), point.len()); + + let cell = geometry.index_map()[cell_index]; + + for component in point.iter_mut() { + *component = T::from(0.0).unwrap(); + } + for (i, v) in geometry.cell_points(cell).unwrap().iter().enumerate() { + let t = unsafe { *table.get_unchecked([0, point_index, i, 0]) }; + for (j, component) in point.iter_mut().enumerate() { + *component += *geometry.coordinate(*v, j).unwrap() * t; + } + } +} + +pub fn compute_jacobian + Shape<4>>( + geometry: &impl Geometry, + table: Table, + tdim: usize, + cell_index: usize, + point_index: usize, + jacobian: &mut [T], +) { + let gdim = geometry.dim(); + assert_eq!(jacobian.len(), gdim * tdim); + + let cell = geometry.index_map()[cell_index]; + + for component in jacobian.iter_mut() { + *component = T::from(0.0).unwrap(); + } + for (i, v) in geometry.cell_points(cell).unwrap().iter().enumerate() { + for gd in 0..gdim { + for td in 0..tdim { + jacobian[td * gdim + gd] += *geometry.coordinate(*v, gd).unwrap() + * unsafe { *table.get_unchecked([1 + td, point_index, i, 0]) }; + } + } + } +} + +pub fn compute_normal_from_jacobian23(jacobian: &[T], normal: &mut [T]) { + assert_eq!(jacobian.len(), 6); + assert_eq!(normal.len(), 3); + + for (i, j, k) in [(0, 1, 2), (1, 2, 0), (2, 0, 1)] { + normal[i] = jacobian[j] * jacobian[3 + k] - jacobian[k] * jacobian[3 + j]; + } + let size = Scalar::sqrt(normal.iter().map(|&x| Scalar::powi(x, 2)).sum::()); + for i in normal.iter_mut() { + *i /= size; + } +} diff --git a/src/grid_impl/flat_triangle_grid.rs b/src/grid_impl/flat_triangle_grid.rs new file mode 100644 index 0000000..af2e95b --- /dev/null +++ b/src/grid_impl/flat_triangle_grid.rs @@ -0,0 +1,5 @@ +mod geometry; +mod grid; +mod topology; + +pub use crate::grid_impl::flat_triangle_grid::grid::SerialFlatTriangleGrid; diff --git a/src/grid_impl/flat_triangle_grid/geometry.rs b/src/grid_impl/flat_triangle_grid/geometry.rs new file mode 100644 index 0000000..827482f --- /dev/null +++ b/src/grid_impl/flat_triangle_grid/geometry.rs @@ -0,0 +1,503 @@ +//! Implementation of grid geometry + +use crate::grid_impl::traits::{Geometry, GeometryEvaluator}; +use crate::reference_cell; +use crate::types::ReferenceCellType; +use bempp_element::element::Inverse; +use bempp_element::element::{create_element, CiarletElement, ElementFamily}; +use bempp_traits::element::{Continuity, FiniteElement}; +use num::Float; +use rlst::rlst_static_array; +use rlst_common::types::Scalar; +use rlst_dense::{ + array::{Array, SliceArray}, + base_array::BaseArray, + data_container::VectorContainer, + rlst_array_from_slice2, + traits::{ + DefaultIterator, DefaultIteratorMut, RandomAccessByRef, RawAccess, Shape, + UnsafeRandomAccessByRef, + }, +}; +use rlst_proc_macro::rlst_static_type; + +/// Geometry of a serial grid +pub struct SerialFlatTriangleGeometry { + dim: usize, + index_map: Vec, + pub(crate) coordinates: Array, 2>, 2>, + pub(crate) cells: Vec>, + pub(crate) element: CiarletElement, + midpoints: Vec, + diameters: Vec, + volumes: Vec, + pub(crate) normals: Vec, + pub(crate) jacobians: Vec, + cell_indices: Vec, +} + +unsafe impl Sync for SerialFlatTriangleGeometry {} + +impl + Inverse> SerialFlatTriangleGeometry { + pub fn new( + coordinates: Array, 2>, 2>, + cells_input: &[usize], + ) -> Self { + let dim = coordinates.shape()[1]; + assert_eq!(dim, 3); + let ncells = cells_input.len() / dim; + + let mut cells = vec![vec![]; ncells]; + for (cell_i, cell) in cells.iter_mut().enumerate() { + cell.extend_from_slice(&cells_input[cell_i * 3..(cell_i + 1) * 3]) + } + + let mut index_map = vec![0; ncells]; + let mut midpoints = vec![]; + let mut diameters = vec![T::from(0.0).unwrap(); ncells]; + let mut volumes = vec![T::from(0.0).unwrap(); ncells]; + let mut normals = vec![]; + let mut jacobians = vec![]; + + let mut a = rlst_static_array!(T, 3); + let mut b = rlst_static_array!(T, 3); + let mut c = rlst_static_array!(T, 3); + + let mut v0 = rlst_static_array!(T, 3); + let mut v1 = rlst_static_array!(T, 3); + let mut v2 = rlst_static_array!(T, 3); + + for (cell_i, index) in index_map.iter_mut().enumerate() { + *index = cell_i; + + midpoints.push(rlst_static_array!(T, 3)); + normals.push(rlst_static_array!(T, 3)); + jacobians.push(rlst_static_array!(T, 3, 2)); + + for (i, c) in v0.iter_mut().enumerate() { + *c = unsafe { *coordinates.get_unchecked([cells[cell_i][0], i]) }; + } + for (i, c) in v1.iter_mut().enumerate() { + *c = unsafe { *coordinates.get_unchecked([cells[cell_i][1], i]) }; + } + for (i, c) in v2.iter_mut().enumerate() { + *c = unsafe { *coordinates.get_unchecked([cells[cell_i][2], i]) }; + } + + midpoints[cell_i].fill_from( + (v0.view() + v1.view() + v2.view()).scalar_mul(T::from(1.0 / 3.0).unwrap()), + ); + + a.fill_from(v1.view() - v0.view()); + b.fill_from(v2.view() - v0.view()); + c.fill_from(v2.view() - v1.view()); + jacobians[cell_i].view_mut().slice(1, 0).fill_from(a.view()); + jacobians[cell_i].view_mut().slice(1, 1).fill_from(b.view()); + + let a_norm = a.view().norm_2(); + let b_norm = b.view().norm_2(); + let c_norm = c.view().norm_2(); + + a.cross(b.view(), normals[cell_i].view_mut()); + + let normal_length = normals[cell_i].view().norm_2(); + normals[cell_i].scale_in_place(T::one() / normal_length); + + volumes[cell_i] = normal_length / T::from(2.0).unwrap(); + + let s = (a_norm + b_norm + c_norm) / T::from(2.0).unwrap(); + + diameters[cell_i] = T::from(2.0).unwrap() + * Scalar::sqrt(((s - a_norm) * (s - b_norm) * (s - c_norm)) / s); + } + + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 1, + Continuity::Continuous, + ); + let cell_indices = (0..ncells).collect::>(); + + Self { + dim, + index_map, + coordinates, + cells, + element, + midpoints, + diameters, + volumes, + normals, + jacobians, + cell_indices, + } + } +} + +impl Geometry for SerialFlatTriangleGeometry { + type IndexType = usize; + type T = T; + type Element = CiarletElement; + type Evaluator<'a> = GeometryEvaluatorFlatTriangle<'a, T> where Self: 'a; + + fn dim(&self) -> usize { + self.dim + } + + fn index_map(&self) -> &[usize] { + &self.index_map + } + + fn coordinate(&self, point_index: usize, coord_index: usize) -> Option<&Self::T> { + self.coordinates.get([point_index, coord_index]) + } + + fn point_count(&self) -> usize { + self.coordinates.shape()[0] + } + + fn cell_points(&self, index: usize) -> Option<&[usize]> { + if index < self.cells.len() { + Some(&self.cells[index]) + } else { + None + } + } + + fn cell_count(&self) -> usize { + self.cells.len() + } + + fn cell_element(&self, index: usize) -> Option<&Self::Element> { + if index < self.cells.len() { + Some(&self.element) + } else { + None + } + } + + fn element_count(&self) -> usize { + 1 + } + fn element(&self, i: usize) -> Option<&Self::Element> { + if i == 0 { + Some(&self.element) + } else { + None + } + } + fn cell_indices(&self, i: usize) -> Option<&[usize]> { + if i == 0 { + Some(&self.cell_indices) + } else { + None + } + } + + fn midpoint(&self, index: usize, point: &mut [Self::T]) { + point.copy_from_slice(self.midpoints[index].data()); + } + + fn diameter(&self, index: usize) -> Self::T { + self.diameters[index] + } + fn volume(&self, index: usize) -> Self::T { + self.volumes[index] + } + + fn get_evaluator<'a>(&'a self, points: &'a [Self::T]) -> GeometryEvaluatorFlatTriangle<'a, T> { + GeometryEvaluatorFlatTriangle::::new(self, points) + } +} + +pub struct GeometryEvaluatorFlatTriangle<'a, T: Float + Scalar + Inverse> { + geometry: &'a SerialFlatTriangleGeometry, + points: SliceArray<'a, T, 2>, +} + +impl<'a, T: Float + Scalar + Inverse> GeometryEvaluatorFlatTriangle<'a, T> { + fn new(geometry: &'a SerialFlatTriangleGeometry, points: &'a [T]) -> Self { + let tdim = reference_cell::dim(geometry.element.cell_type()); + assert_eq!(points.len() % tdim, 0); + let npoints = points.len() / tdim; + Self { + geometry, + points: rlst_array_from_slice2!(T, points, [npoints, tdim]), + } + } +} + +impl<'a, T: Float + Scalar + Inverse> GeometryEvaluator for GeometryEvaluatorFlatTriangle<'a, T> { + type T = T; + + fn point_count(&self) -> usize { + self.points.shape()[0] + } + + fn compute_point(&self, cell_index: usize, point_index: usize, point: &mut [T]) { + let jacobian = &self.geometry.jacobians[cell_index]; + for (index, val_out) in point.iter_mut().enumerate() { + *val_out = self.geometry.coordinates[[self.geometry.cells[cell_index][0], index]] + + jacobian[[index, 0]] * self.points[[point_index, 0]] + + jacobian[[index, 1]] * self.points[[point_index, 1]]; + } + } + + fn compute_jacobian(&self, cell_index: usize, _point_index: usize, jacobian: &mut [T]) { + for (i, j) in jacobian + .iter_mut() + .zip(self.geometry.jacobians[cell_index].iter()) + { + *i = j; + } + } + + fn compute_normal(&self, cell_index: usize, _point_index: usize, normal: &mut [T]) { + for (i, j) in normal + .iter_mut() + .zip(self.geometry.normals[cell_index].iter()) + { + *i = j; + } + } +} + +#[cfg(test)] +mod test { + use crate::grid_impl::flat_triangle_grid::geometry::*; + use approx::*; + use rlst_dense::{ + rlst_dynamic_array2, + traits::{RandomAccessMut, RawAccessMut}, + }; + + fn example_geometry_flat() -> SerialFlatTriangleGeometry { + let mut points = rlst_dynamic_array2!(f64, [4, 3]); + points[[0, 0]] = 0.0; + points[[0, 1]] = 0.0; + points[[0, 2]] = 0.0; + points[[1, 0]] = 1.0; + points[[1, 1]] = 0.0; + points[[1, 2]] = 0.0; + points[[2, 0]] = 1.0; + points[[2, 1]] = 1.0; + points[[2, 2]] = 0.0; + points[[3, 0]] = 0.0; + points[[3, 1]] = 1.0; + points[[3, 2]] = 0.0; + let cells = vec![0, 1, 2, 0, 2, 3]; + SerialFlatTriangleGeometry::new(points, &cells) + } + + fn example_geometry_3d() -> SerialFlatTriangleGeometry { + let mut points = rlst_dynamic_array2!(f64, [4, 3]); + points[[0, 0]] = 0.0; + points[[0, 1]] = 0.0; + points[[0, 2]] = 0.0; + points[[1, 0]] = 1.0; + points[[1, 1]] = 0.0; + points[[1, 2]] = 1.0; + points[[2, 0]] = 1.0; + points[[2, 1]] = 1.0; + points[[2, 2]] = 0.0; + points[[3, 0]] = 0.0; + points[[3, 1]] = 1.0; + points[[3, 2]] = 0.0; + let cells = vec![0, 1, 2, 0, 2, 3]; + SerialFlatTriangleGeometry::new(points, &cells) + } + + #[test] + fn test_counts() { + let g = example_geometry_flat(); + assert_eq!(g.point_count(), 4); + assert_eq!(g.cell_count(), 2); + } + + #[test] + fn test_cell_points() { + let g = example_geometry_flat(); + for (cell_i, points) in [ + vec![ + vec![0.0, 0.0, 0.0], + vec![1.0, 0.0, 0.0], + vec![1.0, 1.0, 0.0], + ], + vec![ + vec![0.0, 0.0, 0.0], + vec![1.0, 1.0, 0.0], + vec![0.0, 1.0, 0.0], + ], + ] + .iter() + .enumerate() + { + let vs = g.cell_points(cell_i).unwrap(); + for (p_i, point) in points.iter().enumerate() { + for (c_i, coord) in point.iter().enumerate() { + assert_relative_eq!( + *coord, + *g.coordinate(vs[p_i], c_i).unwrap(), + epsilon = 1e-12 + ); + } + } + } + } + + fn triangle_points() -> Array, 2>, 2> { + let mut points = rlst_dynamic_array2!(f64, [2, 2]); + *points.get_mut([0, 0]).unwrap() = 0.2; + *points.get_mut([0, 1]).unwrap() = 0.5; + *points.get_mut([1, 0]).unwrap() = 0.6; + *points.get_mut([1, 1]).unwrap() = 0.1; + points + } + + #[test] + fn test_compute_point_flat() { + let g = example_geometry_flat(); + let points = triangle_points(); + + let evaluator = g.get_evaluator(points.data()); + let mut mapped_point = vec![0.0; 3]; + for (cell_i, pts) in [ + vec![vec![0.7, 0.5, 0.0], vec![0.7, 0.1, 0.0]], + vec![vec![0.2, 0.7, 0.0], vec![0.6, 0.7, 0.0]], + ] + .iter() + .enumerate() + { + for (point_i, point) in pts.iter().enumerate() { + evaluator.compute_point(cell_i, point_i, &mut mapped_point); + for (i, j) in mapped_point.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_compute_point_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut mapped_point = vec![0.0; 3]; + for (cell_i, pts) in [ + vec![vec![0.7, 0.5, 0.2], vec![0.7, 0.1, 0.6]], + vec![vec![0.2, 0.7, 0.0], vec![0.6, 0.7, 0.0]], + ] + .iter() + .enumerate() + { + for (point_i, point) in pts.iter().enumerate() { + evaluator.compute_point(cell_i, point_i, &mut mapped_point); + for (i, j) in mapped_point.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_compute_jacobian_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut computed_jacobian = rlst_dynamic_array2!(f64, [3, 2]); + for (cell_i, jacobian) in [ + vec![vec![1.0, 1.0], vec![0.0, 1.0], vec![1.0, 0.0]], + vec![vec![1.0, 0.0], vec![1.0, 1.0], vec![0.0, 0.0]], + ] + .iter() + .enumerate() + { + for point_i in 0..points.shape()[0] { + evaluator.compute_jacobian(cell_i, point_i, computed_jacobian.data_mut()); + for (i, row) in jacobian.iter().enumerate() { + for (j, entry) in row.iter().enumerate() { + assert_relative_eq!( + *entry, + *computed_jacobian.get([i, j]).unwrap(), + epsilon = 1e-12 + ); + } + } + } + } + } + #[test] + fn test_compute_normal_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut computed_normal = vec![0.0; 3]; + for (cell_i, normal) in [ + vec![ + -1.0 / f64::sqrt(3.0), + 1.0 / f64::sqrt(3.0), + 1.0 / f64::sqrt(3.0), + ], + vec![0.0, 0.0, 1.0], + ] + .iter() + .enumerate() + { + for point_i in 0..points.shape()[0] { + evaluator.compute_normal(cell_i, point_i, &mut computed_normal); + assert_relative_eq!( + computed_normal[0] * computed_normal[0] + + computed_normal[1] * computed_normal[1] + + computed_normal[2] * computed_normal[2], + 1.0, + epsilon = 1e-12 + ); + for (i, j) in computed_normal.iter().zip(normal) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_midpoint_flat() { + let g = example_geometry_flat(); + + let mut midpoint = vec![0.0; 3]; + for (cell_i, point) in [ + vec![2.0 / 3.0, 1.0 / 3.0, 0.0], + vec![1.0 / 3.0, 2.0 / 3.0, 0.0], + ] + .iter() + .enumerate() + { + g.midpoint(cell_i, &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + + #[test] + fn test_midpoint_3d() { + let g = example_geometry_3d(); + + let mut midpoint = vec![0.0; 3]; + for (cell_i, point) in [ + vec![2.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], + vec![1.0 / 3.0, 2.0 / 3.0, 0.0], + ] + .iter() + .enumerate() + { + g.midpoint(cell_i, &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } +} diff --git a/src/grid_impl/flat_triangle_grid/grid.rs b/src/grid_impl/flat_triangle_grid/grid.rs new file mode 100644 index 0000000..9aa8f8b --- /dev/null +++ b/src/grid_impl/flat_triangle_grid/grid.rs @@ -0,0 +1,54 @@ +//! Serial implementation of a grid + +use crate::grid_impl::flat_triangle_grid::{ + geometry::SerialFlatTriangleGeometry, topology::SerialFlatTriangleTopology, +}; +use crate::grid_impl::traits::Grid; +use bempp_element::element::Inverse; +use num::Float; +use rlst_common::types::Scalar; +use rlst_dense::{array::Array, base_array::BaseArray, data_container::VectorContainer}; + +/// A serial grid +pub struct SerialFlatTriangleGrid + Inverse> { + topology: SerialFlatTriangleTopology, + geometry: SerialFlatTriangleGeometry, +} + +impl + Inverse> SerialFlatTriangleGrid { + pub fn new(points: Array, 2>, 2>, cells: &[usize]) -> Self { + // Create the topology + let topology = SerialFlatTriangleTopology::new(cells); + + // Create the geometry + let geometry = SerialFlatTriangleGeometry::::new(points, cells); + + Self { topology, geometry } + } +} + +/// A grid +impl + Inverse> Grid for SerialFlatTriangleGrid { + type T = T; + + /// The type that implements [Topology] + type Topology = SerialFlatTriangleTopology; + + /// The type that implements [Geometry] + type Geometry = SerialFlatTriangleGeometry; + + /// Get the grid topology (See [Topology]) + fn topology(&self) -> &Self::Topology { + &self.topology + } + + /// Get the grid geometry (See [Geometry]) + fn geometry(&self) -> &Self::Geometry { + &self.geometry + } + + // Check if the grid is stored in serial + fn is_serial(&self) -> bool { + true + } +} diff --git a/src/grid_impl/flat_triangle_grid/topology.rs b/src/grid_impl/flat_triangle_grid/topology.rs new file mode 100644 index 0000000..9414231 --- /dev/null +++ b/src/grid_impl/flat_triangle_grid/topology.rs @@ -0,0 +1,259 @@ +//! Implementation of grid topology + +use crate::grid_impl::traits::{Ownership, Topology}; +use crate::reference_cell; +use crate::reference_cell::ReferenceCellType; +use crate::types::CellLocalIndexPair; +use std::collections::HashMap; + +/// Topology of a serial grid +pub struct SerialFlatTriangleTopology { + index_map: Vec, + entities_to_vertices: Vec>>, + cells_to_entities: Vec>>, + entities_to_cells: Vec>>>, + entity_types: Vec, +} + +unsafe impl Sync for SerialFlatTriangleTopology {} + +impl SerialFlatTriangleTopology { + pub fn new(cells: &[usize]) -> Self { + let ncells = cells.len() / 3; + let nvertices = *cells.iter().max().unwrap() + 1; + + let mut index_map = vec![0; ncells]; + + let entity_types = vec![ + ReferenceCellType::Point, + ReferenceCellType::Interval, + ReferenceCellType::Triangle, + ]; + + let mut cells_to_entities = vec![vec![vec![]; ncells]; 3]; + let mut entities_to_cells = vec![vec![]; 3]; + let mut entities_to_vertices = vec![vec![]; 2]; + + let mut edge_indices = HashMap::new(); + + entities_to_cells[2] = vec![vec![]; ncells]; + entities_to_vertices[0] = (0..nvertices).map(|i| vec![i]).collect::>(); + entities_to_cells[0] = vec![vec![]; nvertices]; + + for (cell_i, i) in index_map.iter_mut().enumerate() { + let cell = &cells[cell_i * 3..(cell_i + 1) * 3]; + *i = cell_i; + for (local_index, v) in cell.iter().enumerate() { + entities_to_cells[0][*v].push(CellLocalIndexPair::new(cell_i, local_index)); + } + entities_to_cells[2][cell_i] = vec![CellLocalIndexPair::new(cell_i, 0)]; + cells_to_entities[0][cell_i].extend_from_slice(cell); + cells_to_entities[2][cell_i] = vec![cell_i]; + } + + let ref_conn = &reference_cell::connectivity(ReferenceCellType::Triangle)[1]; + for cell_i in 0..ncells { + let cell = &cells[cell_i * 3..(cell_i + 1) * 3]; + for (local_index, rc) in ref_conn.iter().enumerate() { + let mut first = cell[rc[0][0]]; + let mut second = cell[rc[0][1]]; + if first > second { + std::mem::swap(&mut first, &mut second); + } + if let Some(edge_index) = edge_indices.get(&(first, second)) { + cells_to_entities[1][cell_i].push(*edge_index); + entities_to_cells[1][*edge_index] + .push(CellLocalIndexPair::new(cell_i, local_index)); + } else { + edge_indices.insert((first, second), entities_to_vertices[1].len()); + cells_to_entities[1][cell_i].push(entities_to_vertices[1].len()); + entities_to_cells[1].push(vec![CellLocalIndexPair::new(cell_i, local_index)]); + entities_to_vertices[1].push(vec![first, second]); + } + } + } + + Self { + index_map, + entities_to_vertices, + cells_to_entities, + entities_to_cells, + entity_types, + } + } +} + +impl Topology for SerialFlatTriangleTopology { + type IndexType = usize; + + fn dim(&self) -> usize { + 2 + } + fn index_map(&self) -> &[Self::IndexType] { + &self.index_map + } + fn entity_count(&self, etype: ReferenceCellType) -> usize { + if self.entity_types.contains(&etype) { + self.entities_to_cells[reference_cell::dim(etype)].len() + } else { + 0 + } + } + fn entity_count_by_dim(&self, dim: usize) -> usize { + self.entity_count(self.entity_types[dim]) + } + fn cell(&self, index: Self::IndexType) -> Option<&[usize]> { + if index < self.cells_to_entities[2].len() { + Some(&self.cells_to_entities[2][index]) + } else { + None + } + } + fn cell_type(&self, index: Self::IndexType) -> Option { + if index < self.cells_to_entities[2].len() { + Some(self.entity_types[2]) + } else { + None + } + } + + fn entity_types(&self, dim: usize) -> &[ReferenceCellType] { + &self.entity_types[dim..dim + 1] + } + + fn entity_ownership(&self, _dim: usize, _index: Self::IndexType) -> Ownership { + Ownership::Owned + } + fn cell_to_entities(&self, index: Self::IndexType, dim: usize) -> Option<&[Self::IndexType]> { + if dim <= 2 && index < self.cells_to_entities[dim].len() { + Some(&self.cells_to_entities[dim][index]) + } else { + None + } + } + fn entity_to_cells( + &self, + dim: usize, + index: Self::IndexType, + ) -> Option<&[CellLocalIndexPair]> { + if dim <= 2 && index < self.entities_to_cells[dim].len() { + Some(&self.entities_to_cells[dim][index]) + } else { + None + } + } + + fn entity_vertices(&self, dim: usize, index: Self::IndexType) -> Option<&[Self::IndexType]> { + if dim == 2 { + self.cell_to_entities(index, 0) + } else if dim < 2 && index < self.entities_to_vertices[dim].len() { + Some(&self.entities_to_vertices[dim][index]) + } else { + None + } + } +} + +#[cfg(test)] +mod test { + use crate::grid_impl::flat_triangle_grid::topology::*; + + fn example_topology() -> SerialFlatTriangleTopology { + SerialFlatTriangleTopology::new(&[0, 1, 2, 2, 1, 3]) + } + + #[test] + fn test_counts() { + let t = example_topology(); + assert_eq!(t.dim(), 2); + assert_eq!(t.entity_count(ReferenceCellType::Point), 4); + assert_eq!(t.entity_count(ReferenceCellType::Interval), 5); + assert_eq!(t.entity_count(ReferenceCellType::Triangle), 2); + } + + #[test] + fn test_cell_entities_points() { + let t = example_topology(); + for (i, vertices) in [[0, 1, 2], [2, 1, 3]].iter().enumerate() { + let c = t.cell_to_entities(i, 0).unwrap(); + assert_eq!(c.len(), 3); + assert_eq!(c, vertices); + } + } + + #[test] + fn test_cell_entities_intervals() { + let t = example_topology(); + for (i, edges) in [[0, 1, 2], [3, 4, 0]].iter().enumerate() { + let c = t.cell_to_entities(i, 1).unwrap(); + assert_eq!(c.len(), 3); + assert_eq!(c, edges); + } + } + + #[test] + fn test_cell_entities_triangles() { + let t = example_topology(); + for i in 0..2 { + let c = t.cell_to_entities(i, 2).unwrap(); + assert_eq!(c.len(), 1); + assert_eq!(c[0], i); + } + } + + #[test] + fn test_entities_to_cells_points() { + let t = example_topology(); + let c_to_e = (0..t.entity_count(ReferenceCellType::Triangle)) + .map(|i| t.cell_to_entities(i, 0).unwrap()) + .collect::>(); + let e_to_c = (0..t.entity_count(ReferenceCellType::Point)) + .map(|i| { + t.entity_to_cells(0, i) + .unwrap() + .iter() + .map(|x| x.cell) + .collect::>() + }) + .collect::>(); + + for (i, cell) in c_to_e.iter().enumerate() { + for v in *cell { + assert!(e_to_c[*v].contains(&i)); + } + } + for (i, cells) in e_to_c.iter().enumerate() { + for c in cells { + assert!(c_to_e[*c].contains(&i)); + } + } + } + + #[test] + fn test_entities_to_cells_edges() { + let t = example_topology(); + let c_to_e = (0..t.entity_count(ReferenceCellType::Triangle)) + .map(|i| t.cell_to_entities(i, 1).unwrap()) + .collect::>(); + let e_to_c = (0..t.entity_count(ReferenceCellType::Interval)) + .map(|i| { + t.entity_to_cells(1, i) + .unwrap() + .iter() + .map(|x| x.cell) + .collect::>() + }) + .collect::>(); + + for (i, cell) in c_to_e.iter().enumerate() { + for v in *cell { + assert!(e_to_c[*v].contains(&i)); + } + } + for (i, cells) in e_to_c.iter().enumerate() { + for c in cells { + assert!(c_to_e[*c].contains(&i)); + } + } + } +} diff --git a/src/grid_impl/grid.rs b/src/grid_impl/grid.rs index ecde1fc..f683af1 100644 --- a/src/grid_impl/grid.rs +++ b/src/grid_impl/grid.rs @@ -1,12 +1,14 @@ -use crate::grid_impl::traits::{Geometry, Grid, Topology}; +use crate::grid_impl::traits::{Geometry, GeometryEvaluator, Grid, Topology}; +use crate::reference_cell; use crate::reference_cell::ReferenceCellType; use crate::traits::{ cell::CellType, geometry::GeometryType, grid::GridType, point::PointType, reference_map::ReferenceMapType, topology::TopologyType, }; -use crate::types::vertex_iterator::PointIterator; use crate::types::CellLocalIndexPair; +use bempp_traits::element::FiniteElement; use num::Float; +use rlst_common::types::Scalar; use std::iter::Copied; use std::marker::PhantomData; @@ -21,14 +23,40 @@ pub struct Cell<'a, T: Float, GridImpl: Grid> { _t: PhantomData, } pub struct CellTopology<'a, GridImpl: Grid> { - grid: &'a GridImpl, + topology: &'a ::Topology, index: <::Topology as Topology>::IndexType, } pub struct CellGeometry<'a, T: Float, GridImpl: Grid> { - grid: &'a GridImpl, + geometry: &'a ::Geometry, index: <::Geometry as Geometry>::IndexType, _t: PhantomData, } +pub struct ReferenceMap<'a, GridImpl: Grid> { + grid: &'a GridImpl, + evaluator: <::Geometry as Geometry>::Evaluator<'a>, +} +pub struct PointIterator<'a, GridImpl: Grid, Iter: std::iter::Iterator> { + iter: Iter, + geometry: &'a ::Geometry, +} + +impl<'a, GridImpl: Grid, Iter: std::iter::Iterator> std::iter::Iterator + for PointIterator<'a, GridImpl, Iter> +{ + type Item = Point<'a, ::T, ::Geometry>; + + fn next(&mut self) -> Option { + if let Some(index) = self.iter.next() { + Some(Point { + geometry: self.geometry, + index, + _t: PhantomData, + }) + } else { + None + } + } +} impl<'a, T: Float, G: Geometry> PointType for Point<'a, T, G> { type T = T; @@ -46,7 +74,7 @@ impl<'a, T: Float, G: Geometry> PointType for Point<'a, T, G> { } } -impl<'grid, T: Float, GridImpl: Grid> CellType for Cell<'grid, T, GridImpl> +impl<'grid, T: Float + Scalar, GridImpl: Grid> CellType for Cell<'grid, T, GridImpl> where GridImpl: 'grid, { @@ -64,7 +92,7 @@ where fn topology(&self) -> Self::Topology<'_> { CellTopology::<'_, GridImpl> { - grid: self.grid, // TODO: replace with just topology + topology: self.grid.topology(), index: self.grid.topology().index_map()[self.index], } } @@ -75,14 +103,14 @@ where fn geometry(&self) -> Self::Geometry<'_> { CellGeometry::<'_, T, GridImpl> { - grid: self.grid, + geometry: self.grid.geometry(), index: self.grid.geometry().index_map()[self.index], _t: PhantomData, } } } -impl<'grid, GridImpl: Grid> TopologyType for CellTopology<'grid, GridImpl> +impl<'grid, T: Float + Scalar, GridImpl: Grid> TopologyType for CellTopology<'grid, GridImpl> where GridImpl: 'grid, { @@ -99,38 +127,36 @@ where Self: 'a; fn vertex_indices(&self) -> Self::VertexIndexIter<'_> { - self.grid - .topology() - .connectivity(self.grid.topology().dim(), self.index, 0) + self.topology + .cell_to_entities(self.index, 0) .unwrap() .iter() .copied() } fn edge_indices(&self) -> Self::EdgeIndexIter<'_> { - self.grid - .topology() - .connectivity(self.grid.topology().dim(), self.index, 1) + self.topology + .cell_to_entities(self.index, 1) .unwrap() .iter() .copied() } fn face_indices(&self) -> Self::FaceIndexIter<'_> { - self.grid - .topology() - .connectivity(self.grid.topology().dim(), self.index, 2) + self.topology + .cell_to_entities(self.index, 2) .unwrap() .iter() .copied() } fn cell_type(&self) -> ReferenceCellType { - self.grid.topology().cell_type(self.index).unwrap() + self.topology.cell_type(self.index).unwrap() } } -impl<'grid, T: Float, GridImpl: Grid> GeometryType for CellGeometry<'grid, T, GridImpl> +impl<'grid, T: Float + Scalar, GridImpl: Grid> GeometryType + for CellGeometry<'grid, T, GridImpl> where GridImpl: 'grid, { @@ -142,103 +168,90 @@ where type PointIterator<'iter> = Self::VertexIterator<'iter> where Self: 'iter; fn physical_dimension(&self) -> usize { - self.grid.geometry().dim() + self.geometry.dim() } fn midpoint(&self, point: &mut [T]) { - self.grid.geometry().midpoint(self.index, point) + self.geometry.midpoint(self.index, point) } fn diameter(&self) -> T { - self.grid.geometry().diameter(self.index) + self.geometry.diameter(self.index) } fn volume(&self) -> T { - self.grid.geometry().volume(self.index) + self.geometry.volume(self.index) } fn points(&self) -> Self::PointIterator<'_> { - panic!(); + PointIterator { + iter: self + .geometry + .cell_points(self.index) + .unwrap() + .iter() + .copied(), + geometry: self.geometry, + } } fn vertices(&self) -> Self::VertexIterator<'_> { - panic!(); + let cell_type = self.geometry.cell_element(self.index).unwrap().cell_type(); + let nvertices = reference_cell::entity_counts(cell_type)[0]; + PointIterator { + iter: self.geometry.cell_points(self.index).unwrap()[..nvertices] + .iter() + .copied(), + geometry: self.geometry, + } } } -pub struct ReferenceMap<'a, GridImpl: Grid> { - _grid: &'a GridImpl, -} - -impl<'a, GridImpl: Grid> ReferenceMapType for ReferenceMap<'a, GridImpl> { +impl<'a, T: Float + Scalar, GridImpl: Grid> ReferenceMapType for ReferenceMap<'a, GridImpl> { type Grid = GridImpl; fn domain_dimension(&self) -> usize { - panic!(); + self.grid.topology().dim() } fn physical_dimension(&self) -> usize { - panic!(); + self.grid.geometry().dim() } fn number_of_reference_points(&self) -> usize { - panic!(); + self.evaluator.point_count() } fn reference_to_physical( &self, - _point_index: usize, - _value: &mut [::T], + cell_index: usize, + point_index: usize, + value: &mut [::T], ) { - panic!(); + self.evaluator.compute_point(cell_index, point_index, value); } - fn jacobian(&self, _point_index: usize, _value: &mut [::T]) { - panic!(); + fn jacobian(&self, cell_index: usize, point_index: usize, value: &mut [T]) { + self.evaluator + .compute_jacobian(cell_index, point_index, value); } - fn normal(&self, _point_index: usize, _value: &mut [::T]) { - panic!(); + fn normal(&self, cell_index: usize, point_index: usize, value: &mut [T]) { + self.evaluator + .compute_normal(cell_index, point_index, value); } } -pub struct ReferenceMapIterator<'a, Iter: std::iter::Iterator, GridImpl: Grid> -where - GridImpl: 'a, -{ - _grid: &'a GridImpl, - _iter: PhantomData, -} - -impl<'a, Iter: std::iter::Iterator, GridImpl: Grid> Iterator - for ReferenceMapIterator<'a, Iter, GridImpl> -{ - type Item = ReferenceMap<'a, GridImpl>; - - fn next(&mut self) -> Option { - panic!(); - } -} - -impl<'grid, T: Float, GridImpl: Grid> GridType for GridImpl +impl<'grid, T: Float + Scalar, GridImpl: Grid> GridType for GridImpl where GridImpl: 'grid, { type T = T; + type IndexType = <::Topology as Topology>::IndexType; type ReferenceMap<'a> = ReferenceMap<'a, GridImpl> where Self: 'a; - type ReferenceMapIterator<'a, Iter: std::iter::Iterator> = - ReferenceMapIterator<'a, Iter, GridImpl> - where - Self: 'a, - Iter: 'a; - - // PROPOSAL: - // Vertex = one of the corners of a cell - // Point = a point in the geometry - type Point<'a> = Point<'a, T, GridImpl::Geometry> where Self: 'a; type Cell<'a> = Cell<'a, T, GridImpl> where Self: 'a; @@ -284,33 +297,27 @@ where fn reference_to_physical_map<'a>( &'a self, - _reference_points: &'a [Self::T], - _cell_index: usize, + reference_points: &'a [Self::T], ) -> Self::ReferenceMap<'a> { - panic!(); - } - - fn iter_reference_to_physical_map<'a, Iter: std::iter::Iterator + 'a>( - &'a self, - _reference_points: &'a [Self::T], - _iter: Iter, - ) -> Self::ReferenceMapIterator<'a, Iter> - where - Self: 'a, - { - panic!(); + Self::ReferenceMap { + grid: self, + evaluator: self.geometry().get_evaluator(reference_points), + } } - fn point_to_cells(&self, _point_index: usize) -> &[CellLocalIndexPair] { - panic!(); + fn vertex_to_cells( + &self, + vertex_index: Self::IndexType, + ) -> &[CellLocalIndexPair] { + self.topology().entity_to_cells(0, vertex_index).unwrap() } - fn edge_to_cells(&self, _edge_index: usize) -> &[CellLocalIndexPair] { - panic!(); + fn edge_to_cells(&self, edge_index: Self::IndexType) -> &[CellLocalIndexPair] { + self.topology().entity_to_cells(1, edge_index).unwrap() } - fn face_to_cells(&self, _face_index: usize) -> &[CellLocalIndexPair] { - panic!(); + fn face_to_cells(&self, face_index: Self::IndexType) -> &[CellLocalIndexPair] { + self.topology().entity_to_cells(2, face_index).unwrap() } } @@ -319,40 +326,40 @@ mod test { use crate::grid_impl::grid::*; use crate::grid_impl::mixed_grid::SerialMixedGrid; use crate::grid_impl::single_element_grid::SerialSingleElementGrid; + use rlst_dense::{rlst_dynamic_array2, traits::RandomAccessMut}; #[test] fn test_grid_mixed_cell_type() { + let mut points = rlst_dynamic_array2!(f64, [9, 3]); + *points.get_mut([0, 0]).unwrap() = -1.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([0, 2]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = -0.5; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([1, 2]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 0.0; + *points.get_mut([2, 1]).unwrap() = 0.0; + *points.get_mut([2, 2]).unwrap() = 0.0; + *points.get_mut([3, 0]).unwrap() = 1.0; + *points.get_mut([3, 1]).unwrap() = 0.0; + *points.get_mut([3, 2]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 2.0; + *points.get_mut([4, 1]).unwrap() = 0.0; + *points.get_mut([4, 2]).unwrap() = 0.0; + *points.get_mut([5, 0]).unwrap() = -std::f64::consts::FRAC_1_SQRT_2; + *points.get_mut([5, 1]).unwrap() = std::f64::consts::FRAC_1_SQRT_2; + *points.get_mut([5, 2]).unwrap() = 0.0; + *points.get_mut([6, 0]).unwrap() = 0.0; + *points.get_mut([6, 1]).unwrap() = 0.5; + *points.get_mut([6, 2]).unwrap() = 0.0; + *points.get_mut([7, 0]).unwrap() = 0.0; + *points.get_mut([7, 1]).unwrap() = 1.0; + *points.get_mut([7, 2]).unwrap() = 0.0; + *points.get_mut([8, 0]).unwrap() = 1.0; + *points.get_mut([8, 1]).unwrap() = 1.0; + *points.get_mut([8, 2]).unwrap() = 0.0; let grid = SerialMixedGrid::::new( - vec![ - -1.0, - 0.0, - 0.0, - -0.5, - 0.0, - 0.2, - 0.0, - 0.0, - 0.0, - 1.0, - 0.0, - 0.0, - 2.0, - 0.0, - 0.0, - -std::f64::consts::FRAC_1_SQRT_2, - std::f64::consts::FRAC_1_SQRT_2, - 0.0, - 0.0, - 0.5, - 0.0, - 0.0, - 1.0, - 0.0, - 1.0, - 1.0, - 0.0, - ], - 3, + points, &[0, 2, 7, 6, 5, 1, 2, 3, 7, 8, 3, 4, 8], &[ ReferenceCellType::Triangle, @@ -397,12 +404,36 @@ mod test { #[test] fn test_grid_single_element() { + let mut points = rlst_dynamic_array2!(f64, [9, 3]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([0, 2]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.5; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([1, 2]).unwrap() = 0.2; + *points.get_mut([2, 0]).unwrap() = 1.0; + *points.get_mut([2, 1]).unwrap() = 0.0; + *points.get_mut([2, 2]).unwrap() = 0.0; + *points.get_mut([3, 0]).unwrap() = 0.0; + *points.get_mut([3, 1]).unwrap() = 0.5; + *points.get_mut([3, 2]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.5; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([4, 2]).unwrap() = 0.0; + *points.get_mut([5, 0]).unwrap() = 1.0; + *points.get_mut([5, 1]).unwrap() = 0.5; + *points.get_mut([5, 2]).unwrap() = 0.0; + *points.get_mut([6, 0]).unwrap() = 0.0; + *points.get_mut([6, 1]).unwrap() = 1.0; + *points.get_mut([6, 2]).unwrap() = 0.0; + *points.get_mut([7, 0]).unwrap() = 0.5; + *points.get_mut([7, 1]).unwrap() = 1.0; + *points.get_mut([7, 2]).unwrap() = 0.0; + *points.get_mut([8, 0]).unwrap() = 1.0; + *points.get_mut([8, 1]).unwrap() = 1.0; + *points.get_mut([8, 2]).unwrap() = 0.0; let grid = SerialSingleElementGrid::::new( - vec![ - 0.0, 0.0, 0.0, 0.5, 0.0, 0.2, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 1.0, - 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0, - ], - 3, + points, &[0, 2, 6, 4, 3, 1, 2, 8, 6, 7, 4, 5], ReferenceCellType::Triangle, 2, diff --git a/src/grid_impl/mixed_grid/geometry.rs b/src/grid_impl/mixed_grid/geometry.rs index e58efae..39cca48 100644 --- a/src/grid_impl/mixed_grid/geometry.rs +++ b/src/grid_impl/mixed_grid/geometry.rs @@ -1,66 +1,99 @@ //! Implementation of grid geometry -use crate::grid_impl::traits::Geometry; +use crate::grid_impl::common::{compute_jacobian, compute_normal_from_jacobian23, compute_point}; +use crate::grid_impl::traits::{Geometry, GeometryEvaluator}; +use crate::reference_cell; use bempp_element::element::CiarletElement; use bempp_traits::element::FiniteElement; use num::Float; +use rlst_common::types::Scalar; +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_array_from_slice2, rlst_dynamic_array4, + traits::{RandomAccessByRef, Shape, UnsafeRandomAccessByRef}, +}; /// Geometry of a serial grid -pub struct SerialMixedGeometry { +pub struct SerialMixedGeometry { dim: usize, index_map: Vec<(usize, usize)>, - // TODO: change storage to rlst - coordinates: Vec, + coordinates: Array, 2>, 2>, cells: Vec>, - elements: Vec, + pub(crate) elements: Vec>, midpoints: Vec>>, diameters: Vec>, volumes: Vec>, + cell_indices: Vec>, } -unsafe impl Sync for SerialMixedGeometry {} +unsafe impl Sync for SerialMixedGeometry {} -impl SerialMixedGeometry { +impl SerialMixedGeometry { pub fn new( - coordinates: Vec, - dim: usize, + coordinates: Array, 2>, 2>, cells_input: &[usize], - elements: Vec, + elements: Vec>, cell_elements: &[usize], ) -> Self { + let dim = coordinates.shape()[1]; let mut index_map = vec![(0, 0); cell_elements.len()]; let mut cells = vec![]; - let mut midpoints = vec![]; - let mut diameters = vec![]; - let mut volumes = vec![]; + let mut cell_indices = vec![]; for (element_index, _e) in elements.iter().enumerate() { - let mut e_cells = vec![]; let mut start = 0; + cells.push(vec![]); + cell_indices.push(vec![]); + for (cell_i, element_i) in cell_elements.iter().enumerate() { let size = elements[*element_i].dim(); if *element_i == element_index { - index_map[cell_i] = (element_index, e_cells.len() / size); - e_cells.extend_from_slice(&cells_input[start..start + size]); + let cell_index = (element_index, cell_indices[element_index].len()); + index_map[cell_i] = cell_index; + cell_indices[element_index].push(cell_index); + cells[element_index].extend_from_slice(&cells_input[start..start + size]); } start += size; } - cells.push(e_cells); } + let mut midpoints = vec![vec![]; cell_elements.len()]; + let mut diameters = vec![vec![]; cell_elements.len()]; + let mut volumes = vec![vec![]; cell_elements.len()]; + for (element_index, e) in elements.iter().enumerate() { - let mut e_midpoints = vec![]; - let mut e_diameters = vec![]; - let mut e_volumes = vec![]; - for _cell in 0..cells[element_index].len() / e.dim() { - e_midpoints.push(vec![T::from(0.0).unwrap(); dim]); // TODO - e_diameters.push(T::from(0.0).unwrap()); // TODO - e_volumes.push(T::from(0.0).unwrap()); // TODO + let ncells = cells[element_index].len() / e.dim(); + let size = e.dim(); + + midpoints[element_index] = vec![vec![T::from(0.0).unwrap(); dim]; ncells]; + diameters[element_index] = vec![T::from(0.0).unwrap(); ncells]; // TODO + volumes[element_index] = vec![T::from(0.0).unwrap(); ncells]; // TODO + + let mut table = rlst_dynamic_array4!(T, e.tabulate_array_shape(0, 1)); + e.tabulate( + &rlst_array_from_slice2!( + T, + &reference_cell::midpoint(e.cell_type()), + [1, reference_cell::dim(e.cell_type())] + ), + 0, + &mut table, + ); + + let mut start = 0; + for cell_i in 0..ncells { + for (i, v) in cells[element_index][start..start + size].iter().enumerate() { + let t = unsafe { *table.get_unchecked([0, 0, i, 0]) }; + for (j, component) in midpoints[element_index][cell_i].iter_mut().enumerate() { + *component += unsafe { *coordinates.get_unchecked([*v, j]) } * t; + } + } + + start += size; } - midpoints.push(e_midpoints); - diameters.push(e_diameters); - volumes.push(e_volumes); } Self { @@ -72,14 +105,16 @@ impl SerialMixedGeometry { midpoints, diameters, volumes, + cell_indices, } } } -impl Geometry for SerialMixedGeometry { +impl Geometry for SerialMixedGeometry { type IndexType = (usize, usize); type T = T; - type Element = CiarletElement; + type Element = CiarletElement; + type Evaluator<'a> = GeometryEvaluatorMixed<'a, T>; fn dim(&self) -> usize { self.dim @@ -90,15 +125,11 @@ impl Geometry for SerialMixedGeometry { } fn coordinate(&self, point_index: usize, coord_index: usize) -> Option<&Self::T> { - if coord_index < self.dim && point_index * self.dim < self.coordinates.len() { - Some(&self.coordinates[point_index * self.dim + coord_index]) - } else { - None - } + self.coordinates.get([point_index, coord_index]) } fn point_count(&self) -> usize { - self.coordinates.len() / self.dim + self.coordinates.shape()[0] } fn cell_points(&self, index: (usize, usize)) -> Option<&[usize]> { @@ -140,9 +171,9 @@ impl Geometry for SerialMixedGeometry { None } } - fn cells(&self, i: usize) -> Option<&[usize]> { + fn cell_indices(&self, i: usize) -> Option<&[Self::IndexType]> { if i < self.cells.len() { - Some(&self.cells[i]) + Some(&self.cell_indices[i]) } else { None } @@ -158,6 +189,80 @@ impl Geometry for SerialMixedGeometry { fn volume(&self, index: (usize, usize)) -> Self::T { self.volumes[index.0][index.1] } + + fn get_evaluator<'a>(&'a self, points: &'a [Self::T]) -> Self::Evaluator<'a> { + GeometryEvaluatorMixed::new(self, points) + } +} + +pub struct GeometryEvaluatorMixed<'a, T: Float + Scalar> { + geometry: &'a SerialMixedGeometry, + tdim: usize, + tables: Vec, 4>, 4>>, +} + +impl<'a, T: Float + Scalar> GeometryEvaluatorMixed<'a, T> { + fn new(geometry: &'a SerialMixedGeometry, points: &'a [T]) -> Self { + let tdim = reference_cell::dim(geometry.elements[0].cell_type()); + assert_eq!(points.len() % tdim, 0); + let npoints = points.len() / tdim; + let rlst_points = rlst_array_from_slice2!(T, points, [tdim, npoints]); + + let mut tables = vec![]; + for e in &geometry.elements { + assert_eq!(reference_cell::dim(e.cell_type()), tdim); + let mut table = rlst_dynamic_array4!(T, e.tabulate_array_shape(1, npoints)); + e.tabulate(&rlst_points, 1, &mut table); + tables.push(table); + } + Self { + geometry, + tdim, + tables, + } + } +} + +impl<'a, T: Float + Scalar> GeometryEvaluator for GeometryEvaluatorMixed<'a, T> { + type T = T; + + fn point_count(&self) -> usize { + self.tables[0].shape()[1] + } + + fn compute_point(&self, cell_index: usize, point_index: usize, point: &mut [T]) { + let cell = self.geometry.index_map()[cell_index]; + compute_point( + self.geometry, + self.tables[cell.0].view(), + cell_index, + point_index, + point, + ); + } + + fn compute_jacobian(&self, cell_index: usize, point_index: usize, jacobian: &mut [T]) { + let cell = self.geometry.index_map()[cell_index]; + compute_jacobian( + self.geometry, + self.tables[cell.0].view(), + self.tdim, + cell_index, + point_index, + jacobian, + ); + } + + fn compute_normal(&self, cell_index: usize, point_index: usize, normal: &mut [T]) { + let gdim = self.geometry.dim(); + let tdim = self.tdim; + assert_eq!(tdim, 2); + assert_eq!(tdim, gdim - 1); + + let mut jacobian = vec![T::from(0.0).unwrap(); gdim * tdim]; + self.compute_jacobian(cell_index, point_index, &mut jacobian[..]); + compute_normal_from_jacobian23(&jacobian, normal); + } } #[cfg(test)] @@ -166,6 +271,7 @@ mod test { use approx::*; use bempp_element::element::{create_element, ElementFamily}; use bempp_traits::element::Continuity; + use rlst_dense::{rlst_dynamic_array2, traits::RandomAccessMut}; fn example_geometry() -> SerialMixedGeometry { let p1triangle = create_element( @@ -174,12 +280,47 @@ mod test { 1, Continuity::Continuous, ); + let mut points = rlst_dynamic_array2!(f64, [4, 2]); + *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, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 1.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + SerialMixedGeometry::new(points, &[0, 1, 2, 0, 2, 3], vec![p1triangle], &[0, 0]) + } + + fn example_geometry_mixed() -> SerialMixedGeometry { + let p1triangle = create_element( + ElementFamily::Lagrange, + bempp_element::cell::ReferenceCellType::Triangle, + 1, + Continuity::Continuous, + ); + let p1quad = create_element( + ElementFamily::Lagrange, + bempp_element::cell::ReferenceCellType::Quadrilateral, + 1, + Continuity::Continuous, + ); + let mut points = rlst_dynamic_array2!(f64, [5, 2]); + *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, 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() = 2.0; + *points.get_mut([4, 1]).unwrap() = 0.0; SerialMixedGeometry::new( - vec![0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0], - 2, - &[0, 1, 2, 0, 2, 3], - vec![p1triangle], - &[0, 0], + points, + &[0, 1, 2, 3, 1, 4, 3], + vec![p1quad, p1triangle], + &[0, 1], ) } @@ -203,34 +344,16 @@ mod test { let vs = g.cell_points((0, cell_i)).unwrap(); for (p_i, point) in points.iter().enumerate() { for (c_i, coord) in point.iter().enumerate() { - assert_relative_eq!(*coord, *g.coordinate(vs[p_i], c_i).unwrap()); + assert_relative_eq!( + *coord, + *g.coordinate(vs[p_i], c_i).unwrap(), + epsilon = 1e-12 + ); } } } } - fn example_geometry_mixed() -> SerialMixedGeometry { - let p1triangle = create_element( - ElementFamily::Lagrange, - bempp_element::cell::ReferenceCellType::Triangle, - 1, - Continuity::Continuous, - ); - let p1quad = create_element( - ElementFamily::Lagrange, - bempp_element::cell::ReferenceCellType::Quadrilateral, - 1, - Continuity::Continuous, - ); - SerialMixedGeometry::new( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 0.0], - 2, - &[0, 1, 2, 3, 1, 4, 3], - vec![p1quad, p1triangle], - &[0, 1], - ) - } - #[test] fn test_counts_mixed() { let g = example_geometry_mixed(); @@ -256,9 +379,45 @@ mod test { let vs = g.cell_points(cell_i).unwrap(); for (p_i, point) in points.iter().enumerate() { for (c_i, coord) in point.iter().enumerate() { - assert_relative_eq!(*coord, *g.coordinate(vs[p_i], c_i).unwrap()); + assert_relative_eq!( + *coord, + *g.coordinate(vs[p_i], c_i).unwrap(), + epsilon = 1e-12 + ); } } } } + + #[test] + fn test_midpoint() { + let g = example_geometry(); + + let mut midpoint = vec![0.0; 2]; + for (cell_i, point) in [vec![2.0 / 3.0, 1.0 / 3.0], vec![1.0 / 3.0, 2.0 / 3.0]] + .iter() + .enumerate() + { + g.midpoint(g.index_map()[cell_i], &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + + #[test] + fn test_midpoint_mixed() { + let g = example_geometry_mixed(); + + let mut midpoint = vec![0.0; 2]; + for (cell_i, point) in [vec![0.5, 0.5], vec![4.0 / 3.0, 1.0 / 3.0]] + .iter() + .enumerate() + { + g.midpoint(g.index_map()[cell_i], &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } } diff --git a/src/grid_impl/mixed_grid/grid.rs b/src/grid_impl/mixed_grid/grid.rs index c476daa..1c1dcdb 100644 --- a/src/grid_impl/mixed_grid/grid.rs +++ b/src/grid_impl/mixed_grid/grid.rs @@ -4,20 +4,22 @@ use crate::grid_impl::mixed_grid::{geometry::SerialMixedGeometry, topology::Seri use crate::grid_impl::traits::Grid; use crate::reference_cell; use crate::reference_cell::ReferenceCellType; -use bempp_element::element::{create_element, ElementFamily}; +use bempp_element::element::{create_element, ElementFamily, Inverse}; use bempp_traits::element::{Continuity, FiniteElement}; +use log::warn; use num::Float; +use rlst_common::types::Scalar; +use rlst_dense::{array::Array, base_array::BaseArray, data_container::VectorContainer}; /// A serial grid -pub struct SerialMixedGrid { +pub struct SerialMixedGrid { topology: SerialMixedTopology, geometry: SerialMixedGeometry, } -impl SerialMixedGrid { +impl SerialMixedGrid { pub fn new( - points: Vec, - gdim: usize, + points: Array, 2>, 2>, cells: &[usize], cell_types: &[ReferenceCellType], cell_degrees: &[usize], @@ -36,32 +38,12 @@ impl SerialMixedGrid { let elements = element_info .iter() .map(|(i, j)| { - create_element( - ElementFamily::Lagrange, - // TODO: remove this match once bempp-rs and grid-rs use the same ReferenceCellType - match i { - ReferenceCellType::Interval => { - bempp_element::cell::ReferenceCellType::Interval - } - ReferenceCellType::Triangle => { - bempp_element::cell::ReferenceCellType::Triangle - } - ReferenceCellType::Quadrilateral => { - bempp_element::cell::ReferenceCellType::Quadrilateral - } - _ => { - panic!("Unsupported cell type: {:?}", i); - } - }, - *j, - Continuity::Continuous, - ) + create_element::(ElementFamily::Lagrange, *i, *j, Continuity::Continuous) }) .collect::>(); if elements.len() == 1 { - println!("Creating a mixed grid with only one element. Using a SerialSingleElementGrid would be faster."); - // TODO: make into a warning + warn!("Creating a mixed grid with only one element. Using a SerialSingleElementGrid would be faster."); } let mut cell_vertices = vec![]; @@ -78,15 +60,14 @@ impl SerialMixedGrid { let topology = SerialMixedTopology::new(&cell_vertices, cell_types); // Create the geometry - let geometry = - SerialMixedGeometry::::new(points, gdim, cells, elements, &element_numbers); + let geometry = SerialMixedGeometry::::new(points, cells, elements, &element_numbers); Self { topology, geometry } } } /// A grid -impl Grid for SerialMixedGrid { +impl Grid for SerialMixedGrid { type T = T; /// The type that implements [Topology] diff --git a/src/grid_impl/mixed_grid/topology.rs b/src/grid_impl/mixed_grid/topology.rs index 22cc67d..99bc7d2 100644 --- a/src/grid_impl/mixed_grid/topology.rs +++ b/src/grid_impl/mixed_grid/topology.rs @@ -3,12 +3,10 @@ use crate::grid_impl::traits::{Ownership, Topology}; use crate::reference_cell; use crate::reference_cell::ReferenceCellType; +use crate::types::CellLocalIndexPair; use std::collections::HashMap; -// TODO: use 2D rlst arrays here where possible -type Connectivity = Vec>; - fn all_equal(a: &[T], b: &[T]) -> bool { if a.len() != b.len() { false @@ -26,13 +24,15 @@ fn all_in(a: &[T], b: &[T]) -> bool { true } +type IndexType = (ReferenceCellType, usize); + /// Topology of a serial grid pub struct SerialMixedTopology { dim: usize, - index_map: Vec<(ReferenceCellType, usize)>, - cells: HashMap>>, // TODO: use 2D array - connectivity: HashMap>, - cell_connectivity: HashMap, + index_map: Vec, + entities_to_vertices: Vec>>>, + cells_to_entities: Vec>>>, + entities_to_cells: Vec>>>>, entity_types: Vec>, } @@ -45,9 +45,10 @@ impl SerialMixedTopology { let dim = reference_cell::dim(cell_types[0]); let mut entity_types = vec![vec![]; 4]; - let mut cells = HashMap::new(); - let mut connectivity = HashMap::new(); - let mut cell_connectivity = HashMap::new(); + + let mut entities_to_vertices = vec![HashMap::new(); dim]; + let mut cells_to_entities = vec![HashMap::new(); dim + 1]; + let mut entities_to_cells = vec![HashMap::new(); dim + 1]; for c in cell_types { if dim != reference_cell::dim(*c) { @@ -57,7 +58,15 @@ impl SerialMixedTopology { for e in etypes { if !entity_types[dim0].contains(e) { entity_types[dim0].push(*e); - connectivity.insert(*e, vec![vec![]; dim + 1]); + + entities_to_cells[dim0].insert(*e, vec![]); + if dim0 == dim { + for ce in cells_to_entities.iter_mut() { + ce.insert(*e, vec![]); + } + } else { + entities_to_vertices[dim0].insert(*e, vec![]); + } } } } @@ -65,18 +74,19 @@ impl SerialMixedTopology { // dim0 = dim, dim1 = 0 for c in &entity_types[dim] { - let mut subcells = vec![]; - let mut cty = vec![]; - let mut cell_cty = vec![vec![]; dim + 1]; let n = reference_cell::entity_counts(*c)[0]; let mut start = 0; for (i, ct) in cell_types.iter().enumerate() { if *ct == *c { let cell = &cells_input[start..start + n]; - index_map[i] = (*c, subcells.len()); + let cell_i = (*c, cells_to_entities[0][c].len()); + index_map[i] = cell_i; let mut row = vec![]; for v in cell { if !vertices.contains(v) { + for (_, ec) in entities_to_cells[0].iter_mut() { + ec.push(vec![]); + } vertices.push(*v); } row.push(( @@ -84,144 +94,72 @@ impl SerialMixedTopology { vertices.iter().position(|&r| r == *v).unwrap(), )); } - cell_cty[0].extend_from_slice(&row); - subcells.push(row.iter().map(|x| (*c, x.1)).collect()); - cty.push(row); + + for (local_index, v) in row.iter().enumerate() { + entities_to_cells[0].get_mut(&v.0).unwrap()[v.1] + .push(CellLocalIndexPair::new(cell_i, local_index)); + } + entities_to_cells[dim] + .get_mut(c) + .unwrap() + .push(vec![CellLocalIndexPair::new(cell_i, 0)]); + + cells_to_entities[0].get_mut(c).unwrap().push(row); + cells_to_entities[dim] + .get_mut(c) + .unwrap() + .push(vec![cell_i]); } start += reference_cell::entity_counts(*ct)[0]; } - connectivity.get_mut(c).unwrap()[0] = cty; - cell_connectivity.insert(*c, cell_cty); - cells.insert(*c, subcells); + } + for i in 0..vertices.len() { + entities_to_vertices[0] + .get_mut(&ReferenceCellType::Point) + .unwrap() + .push(vec![(ReferenceCellType::Point, i)]); } - // dim1 == 0 - for (dim0, etypes0) in entity_types.iter().enumerate().take(dim) { + for (d, etypes0) in entity_types.iter().enumerate().take(dim).skip(1) { for etype in etypes0 { - let mut cty: Vec> = vec![]; for cell_type in &entity_types[dim] { - let ref_conn = &reference_cell::connectivity(*cell_type)[dim0]; - for cell in &connectivity[&cell_type][0] { - for rc in ref_conn { + let mut c_to_e = vec![]; + let ref_conn = &reference_cell::connectivity(*cell_type)[d]; + for (cell_i, cell) in cells_to_entities[0][&cell_type].iter().enumerate() { + let mut entity_ids = vec![]; + + for (local_index, rc) in ref_conn.iter().enumerate() { let vertices = rc[0].iter().map(|x| cell[*x]).collect::>(); let mut found = false; - for entity in &cty { + for (entity_index, entity) in + entities_to_vertices[d][etype].iter().enumerate() + { if all_equal(entity, &vertices) { + entity_ids.push((*etype, entity_index)); + entities_to_cells[d].get_mut(etype).unwrap()[entity_index] + .push(CellLocalIndexPair::new( + (*cell_type, cell_i), + local_index, + )); + found = true; break; } } if !found { - cty.push(vertices); - } - } - } - } - connectivity.get_mut(etype).unwrap()[0] = cty; - } - } - - // dim0 == dim1 == 0 - let mut nvertices = 0; - let mut cty = vec![]; - for cell_type in &entity_types[dim] { - for cell in &connectivity[&cell_type][0] { - nvertices = std::cmp::max(nvertices, 1 + cell.iter().map(|j| j.1).max().unwrap()); - } - } - for i in 0..nvertices { - cty.push(vec![(ReferenceCellType::Point, i)]); - } - connectivity.get_mut(&ReferenceCellType::Point).unwrap()[0] = cty; - - // dim0 == dim1 - for (dim0, etypes) in entity_types.iter().enumerate().skip(1) { - for etype in etypes { - let mut cty = vec![]; - for i in 0..connectivity[etype][0].len() { - cty.push(vec![(*etype, i)]); - } - connectivity.get_mut(etype).unwrap()[dim0] = cty; - } - } - // dim0 == dim1 == dim - for etype in &entity_types[dim] { - let mut cty = vec![]; - for i in 0..connectivity[etype][0].len() { - cty.push((*etype, i)); - } - cell_connectivity.get_mut(etype).unwrap()[dim] = cty; - } - // dim0 == dim - for cell_type in &entity_types[dim] { - for (dim1, etypes0) in entity_types.iter().enumerate().take(dim).skip(1) { - let mut cty = vec![]; - let mut cell_cty = vec![]; - let entities0 = &connectivity[cell_type][0]; - let ref_conn = &reference_cell::connectivity(*cell_type)[dim1]; - for etype in etypes0 { - let entities1 = &connectivity[etype][0]; - - for entity0 in entities0 { - let mut row = vec![]; - for rc in ref_conn { - let vertices = rc[0].iter().map(|x| entity0[*x]).collect::>(); - for (j, entity1) in entities1.iter().enumerate() { - if all_equal(&vertices, entity1) { - row.push((*etype, j)); - break; - } - } - } - cell_cty.extend_from_slice(&row); - cty.push(row); - } - } - cell_connectivity.get_mut(cell_type).unwrap()[dim1] = cell_cty; - connectivity.get_mut(cell_type).unwrap()[dim1] = cty; - } - } - // dim1 < dim0 - for (dim0, etypes0) in entity_types.iter().enumerate().take(dim + 1).skip(2) { - for etype0 in etypes0 { - for (dim1, etypes1) in entity_types.iter().enumerate().take(dim0).skip(1) { - let mut cty = vec![]; - let entities0 = &connectivity[etype0][0]; - let ref_conn = &reference_cell::connectivity(*etype0)[dim1]; - for etype1 in etypes1 { - let entities1 = &connectivity[etype1][0]; - for entity0 in entities0 { - let mut row = vec![]; - for rc in ref_conn { - let vertices = - rc[0].iter().map(|x| entity0[*x]).collect::>(); - for (j, entity1) in entities1.iter().enumerate() { - if all_equal(&vertices, entity1) { - row.push((*etype1, j)); - break; - } - } + entity_ids.push((*etype, entities_to_vertices[d][etype].len())); + entities_to_cells[d].get_mut(etype).unwrap().push(vec![ + CellLocalIndexPair::new((*cell_type, cell_i), local_index), + ]); + entities_to_vertices[d] + .get_mut(etype) + .unwrap() + .push(vertices); } - cty.push(row); } + c_to_e.push(entity_ids); } - connectivity.get_mut(etype0).unwrap()[dim1] = cty; - } - } - } - // dim1 > dim0 - for (dim0, etypes0) in entity_types.iter().enumerate().take(dim) { - for etype0 in etypes0 { - for (dim1, etypes1) in entity_types.iter().enumerate().take(dim + 1).skip(1) { - let mut data = vec![vec![]; connectivity[etype0][0].len()]; - for etype1 in etypes1 { - for (i, row) in connectivity[etype1][dim0].iter().enumerate() { - for v in row { - data[v.1].push((*etype1, i)); - } - } - } - connectivity.get_mut(etype0).unwrap()[dim1] = data; + *cells_to_entities[d].get_mut(cell_type).unwrap() = c_to_e; } } } @@ -229,16 +167,16 @@ impl SerialMixedTopology { Self { dim, index_map, - cells, - connectivity, - cell_connectivity, + entities_to_vertices, + cells_to_entities, + entities_to_cells, entity_types, } } } impl Topology for SerialMixedTopology { - type IndexType = (ReferenceCellType, usize); + type IndexType = IndexType; fn dim(&self) -> usize { self.dim @@ -247,7 +185,12 @@ impl Topology for SerialMixedTopology { &self.index_map } fn entity_count(&self, etype: ReferenceCellType) -> usize { - self.connectivity[&etype][0].len() + let dim = reference_cell::dim(etype); + if self.entity_types[dim].contains(&etype) { + self.entities_to_cells[dim][&etype].len() + } else { + 0 + } } fn entity_count_by_dim(&self, dim: usize) -> usize { self.entity_types[dim] @@ -255,15 +198,19 @@ impl Topology for SerialMixedTopology { .map(|e| self.entity_count(*e)) .sum() } - fn cell(&self, index: Self::IndexType) -> Option<&[(ReferenceCellType, usize)]> { - if self.cells.contains_key(&index.0) && index.1 < self.cells[&index.0].len() { - Some(&self.cells[&index.0][index.1]) + fn cell(&self, index: Self::IndexType) -> Option<&[IndexType]> { + if self.cells_to_entities[0].contains_key(&index.0) + && index.1 < self.cells_to_entities[0][&index.0].len() + { + Some(&self.cells_to_entities[0][&index.0][index.1]) } else { None } } fn cell_type(&self, index: Self::IndexType) -> Option { - if self.cells.contains_key(&index.0) && index.1 < self.cells[&index.0].len() { + if self.cells_to_entities[0].contains_key(&index.0) + && index.1 < self.cells_to_entities[0][&index.0].len() + { Some(index.0) } else { None @@ -274,46 +221,45 @@ impl Topology for SerialMixedTopology { &self.entity_types[dim] } - fn cell_entities( + fn entity_ownership(&self, _dim: usize, _index: Self::IndexType) -> Ownership { + Ownership::Owned + } + + fn entity_to_cells( &self, - cell_type: ReferenceCellType, dim: usize, - ) -> Option<&[(ReferenceCellType, usize)]> { - if self.cell_connectivity.contains_key(&cell_type) - && dim < self.cell_connectivity[&cell_type].len() + index: Self::IndexType, + ) -> Option<&[CellLocalIndexPair]> { + if dim <= self.dim + && self.entities_to_cells[dim].contains_key(&index.0) + && index.1 < self.entities_to_cells[dim][&index.0].len() { - Some(&self.cell_connectivity[&cell_type][dim]) + Some(&self.entities_to_cells[dim][&index.0][index.1]) } else { None } } - - fn connectivity( - &self, - _dim0: usize, - index: (ReferenceCellType, usize), - dim1: usize, - ) -> Option<&[Self::IndexType]> { - if self.connectivity.contains_key(&index.0) - && dim1 < self.connectivity[&index.0].len() - && index.1 < self.connectivity[&index.0][dim1].len() + fn cell_to_entities(&self, index: Self::IndexType, dim: usize) -> Option<&[Self::IndexType]> { + if dim <= self.dim + && self.cells_to_entities[dim].contains_key(&index.0) + && index.1 < self.cells_to_entities[dim][&index.0].len() { - Some(&self.connectivity[&index.0][dim1][index.1]) + Some(&self.cells_to_entities[dim][&index.0][index.1]) } else { None } } - fn entity_ownership(&self, _dim: usize, _index: Self::IndexType) -> Ownership { - Ownership::Owned - } - - fn extract_index(&self, index: Self::IndexType) -> usize { - let dim = reference_cell::dim(index.0); - if dim < 2 { - index.1 + fn entity_vertices(&self, dim: usize, index: Self::IndexType) -> Option<&[Self::IndexType]> { + if dim == self.dim { + self.cell_to_entities(index, 0) + } else if dim < self.dim + && self.entities_to_vertices[dim].contains_key(&index.0) + && index.1 < self.entities_to_vertices[dim][&index.0].len() + { + Some(&self.entities_to_vertices[dim][&index.0][index.1]) } else { - panic!(); + None } } } @@ -336,133 +282,47 @@ mod test { } #[test] - fn test_cell_entities_points() { - let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 0).unwrap(); - assert_eq!(c.len(), 6); - for i in c { - assert_eq!(i.0, ReferenceCellType::Point); - } - // cell 0 - assert_eq!(c[0].1, 0); - assert_eq!(c[1].1, 1); - assert_eq!(c[2].1, 2); - // cell 1 - assert_eq!(c[3].1, 2); - assert_eq!(c[4].1, 1); - assert_eq!(c[5].1, 3); - } - - #[test] - fn test_cell_entities_intervals() { - let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 1).unwrap(); - assert_eq!(c.len(), 6); - for i in c { - assert_eq!(i.0, ReferenceCellType::Interval); - } - // cell 0 - assert_eq!(c[0].1, 0); - assert_eq!(c[1].1, 1); - assert_eq!(c[2].1, 2); - // cell 1 - assert_eq!(c[3].1, 3); - assert_eq!(c[4].1, 4); - assert_eq!(c[5].1, 0); - } - #[test] - fn test_cell_entities_triangles() { - let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 2).unwrap(); - assert_eq!(c.len(), 2); - // cell 0 - assert_eq!(c[0].0, ReferenceCellType::Triangle); - assert_eq!(c[0].1, 0); - // cell 1 - assert_eq!(c[1].0, ReferenceCellType::Triangle); - assert_eq!(c[1].1, 1); - } - #[test] - fn test_vertex_connectivity() { + fn test_cell_to_entities_vertices() { let t = example_topology(); - - for (id, vertices) in [vec![0], vec![1], vec![2], vec![3]].iter().enumerate() { + for (i, vertices) in [[0, 1, 2], [2, 1, 3]].iter().enumerate() { let c = t - .connectivity(0, (ReferenceCellType::Point, id), 0) + .cell_to_entities((ReferenceCellType::Triangle, i), 0) .unwrap(); - for (i, j) in c.iter().zip(vertices) { + assert_eq!(c.len(), 3); + for i in c { assert_eq!(i.0, ReferenceCellType::Point); - assert_eq!(i.1, *j); - } - } - - for (id, edges) in [vec![1, 2], vec![0, 2, 3], vec![0, 1, 4], vec![3, 4]] - .iter() - .enumerate() - { - let c = t - .connectivity(0, (ReferenceCellType::Point, id), 1) - .unwrap(); - for (i, j) in c.iter().zip(edges) { - assert_eq!(i.0, ReferenceCellType::Interval); - assert_eq!(i.1, *j); - } - } - - for (id, faces) in [vec![0], vec![0, 1], vec![0, 1], vec![1]] - .iter() - .enumerate() - { - let c = t - .connectivity(0, (ReferenceCellType::Point, id), 2) - .unwrap(); - for (i, j) in c.iter().zip(faces) { - assert_eq!(i.0, ReferenceCellType::Triangle); - assert_eq!(i.1, *j); } + assert_eq!(c[0].1, vertices[0]); + assert_eq!(c[1].1, vertices[1]); + assert_eq!(c[2].1, vertices[2]); } } #[test] - fn test_edge_connectivity() { + fn test_cell_to_entities_intervals() { let t = example_topology(); - - for (id, vertices) in [vec![1, 2], vec![0, 2], vec![0, 1], vec![1, 3], vec![2, 3]] - .iter() - .enumerate() - { - let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 0) - .unwrap(); - for (i, j) in c.iter().zip(vertices) { - assert_eq!(i.0, ReferenceCellType::Point); - assert_eq!(i.1, *j); - } - } - - for (id, edges) in [vec![0], vec![1], vec![2], vec![3], vec![4]] - .iter() - .enumerate() - { + for (i, edges) in [[0, 1, 2], [3, 4, 0]].iter().enumerate() { let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 1) + .cell_to_entities((ReferenceCellType::Triangle, i), 1) .unwrap(); - for (i, j) in c.iter().zip(edges) { + assert_eq!(c.len(), 3); + for i in c { assert_eq!(i.0, ReferenceCellType::Interval); - assert_eq!(i.1, *j); } + assert_eq!(c[0].1, edges[0]); + assert_eq!(c[1].1, edges[1]); + assert_eq!(c[2].1, edges[2]); } - - for (id, faces) in [vec![0, 1], vec![0], vec![0], vec![1], vec![1]] - .iter() - .enumerate() - { + } + #[test] + fn test_cell_to_entities_triangles() { + let t = example_topology(); + for i in 0..2 { let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 2) + .cell_to_entities((ReferenceCellType::Triangle, i), 2) .unwrap(); - for (i, j) in c.iter().zip(faces) { - assert_eq!(i.0, ReferenceCellType::Triangle); - assert_eq!(i.1, *j); - } + assert_eq!(c.len(), 1); + assert_eq!(c[0].0, ReferenceCellType::Triangle); + assert_eq!(c[0].1, i); } } @@ -490,7 +350,7 @@ mod test { fn test_mixed_cell_entities_points() { let t = example_topology_mixed(); let c = t - .cell_entities(ReferenceCellType::Quadrilateral, 0) + .cell_to_entities((ReferenceCellType::Quadrilateral, 0), 0) .unwrap(); assert_eq!(c.len(), 4); for i in c { @@ -502,7 +362,9 @@ mod test { assert_eq!(c[2].1, 2); assert_eq!(c[3].1, 3); - let c = t.cell_entities(ReferenceCellType::Triangle, 0).unwrap(); + let c = t + .cell_to_entities((ReferenceCellType::Triangle, 0), 0) + .unwrap(); assert_eq!(c.len(), 3); // cell 1 assert_eq!(c[0].1, 1); @@ -514,7 +376,7 @@ mod test { fn test_mixed_cell_entities_intervals() { let t = example_topology_mixed(); let c = t - .cell_entities(ReferenceCellType::Quadrilateral, 1) + .cell_to_entities((ReferenceCellType::Quadrilateral, 0), 1) .unwrap(); assert_eq!(c.len(), 4); @@ -527,7 +389,9 @@ mod test { assert_eq!(c[2].1, 2); assert_eq!(c[3].1, 3); - let c = t.cell_entities(ReferenceCellType::Triangle, 1).unwrap(); + let c = t + .cell_to_entities((ReferenceCellType::Triangle, 0), 1) + .unwrap(); assert_eq!(c.len(), 3); // cell 1 assert_eq!(c[0].1, 4); @@ -538,158 +402,19 @@ mod test { fn test_mixed_cell_entities_triangles() { let t = example_topology_mixed(); let c = t - .cell_entities(ReferenceCellType::Quadrilateral, 2) + .cell_to_entities((ReferenceCellType::Quadrilateral, 0), 2) .unwrap(); assert_eq!(c.len(), 1); // cell 0 assert_eq!(c[0].0, ReferenceCellType::Quadrilateral); assert_eq!(c[0].1, 0); - let c = t.cell_entities(ReferenceCellType::Triangle, 2).unwrap(); + let c = t + .cell_to_entities((ReferenceCellType::Triangle, 0), 2) + .unwrap(); assert_eq!(c.len(), 1); // cell 1 assert_eq!(c[0].0, ReferenceCellType::Triangle); assert_eq!(c[0].1, 0); } - #[test] - fn test_mixed_vertex_connectivity() { - let t = example_topology_mixed(); - - for (id, vertices) in [vec![0], vec![1], vec![2], vec![3], vec![4]] - .iter() - .enumerate() - { - let c = t - .connectivity(0, (ReferenceCellType::Point, id), 0) - .unwrap(); - for (i, j) in c.iter().zip(vertices) { - assert_eq!(i.0, ReferenceCellType::Point); - assert_eq!(i.1, *j); - } - } - - for (id, edges) in [ - vec![0, 1], - vec![0, 2, 5], - vec![1, 3], - vec![2, 3, 4], - vec![4, 5], - ] - .iter() - .enumerate() - { - let c = t - .connectivity(0, (ReferenceCellType::Point, id), 1) - .unwrap(); - for (i, j) in c.iter().zip(edges) { - assert_eq!(i.0, ReferenceCellType::Interval); - assert_eq!(i.1, *j); - } - } - - for (id, faces) in [ - vec![(ReferenceCellType::Quadrilateral, 0)], - vec![ - (ReferenceCellType::Quadrilateral, 0), - (ReferenceCellType::Triangle, 0), - ], - vec![(ReferenceCellType::Quadrilateral, 0)], - vec![ - (ReferenceCellType::Quadrilateral, 0), - (ReferenceCellType::Triangle, 0), - ], - vec![(ReferenceCellType::Triangle, 0)], - ] - .iter() - .enumerate() - { - let c = t - .connectivity(0, (ReferenceCellType::Point, id), 2) - .unwrap(); - assert_eq!(c, faces); - } - } - #[test] - fn test_mixed_edge_connectivity() { - let t = example_topology_mixed(); - - for (id, vertices) in [ - vec![0, 1], - vec![0, 2], - vec![1, 3], - vec![2, 3], - vec![4, 3], - vec![1, 4], - ] - .iter() - .enumerate() - { - let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 0) - .unwrap(); - for (i, j) in c.iter().zip(vertices) { - assert_eq!(i.0, ReferenceCellType::Point); - assert_eq!(i.1, *j); - } - } - - for (id, edges) in [vec![0], vec![1], vec![2], vec![3], vec![4], vec![5]] - .iter() - .enumerate() - { - let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 1) - .unwrap(); - for (i, j) in c.iter().zip(edges) { - assert_eq!(i.0, ReferenceCellType::Interval); - assert_eq!(i.1, *j); - } - } - - for (id, faces) in [ - vec![(ReferenceCellType::Quadrilateral, 0)], - vec![(ReferenceCellType::Quadrilateral, 0)], - vec![ - (ReferenceCellType::Quadrilateral, 0), - (ReferenceCellType::Triangle, 0), - ], - vec![(ReferenceCellType::Quadrilateral, 0)], - vec![(ReferenceCellType::Triangle, 0)], - vec![(ReferenceCellType::Triangle, 0)], - ] - .iter() - .enumerate() - { - let c = t - .connectivity(1, (ReferenceCellType::Interval, id), 2) - .unwrap(); - assert_eq!(c, faces); - } - } - - fn cell_entities_vs_connectivity(t: &SerialMixedTopology) { - for cell_type in t.entity_types(t.dim()) { - for dim in 0..t.dim() + 1 { - let ce = t.cell_entities(*cell_type, dim).unwrap(); - let n = reference_cell::entity_counts(*cell_type)[dim]; - for i in 0..ce.len() / n { - // TODO: entity_count instead? - let con = t.connectivity(t.dim(), (*cell_type, i), dim).unwrap(); - assert_eq!(con, &ce[n * i..n * (i + 1)]); - } - } - } - } - - #[test] - fn test_cell_entities_vs_connectivity() { - let t = example_topology(); - cell_entities_vs_connectivity(&t); - } - - #[test] - fn test_cell_entities_vs_connectivity_mixed() { - let t = example_topology_mixed(); - cell_entities_vs_connectivity(&t); - } } diff --git a/src/grid_impl/single_element_grid/geometry.rs b/src/grid_impl/single_element_grid/geometry.rs index ea66183..63f3699 100644 --- a/src/grid_impl/single_element_grid/geometry.rs +++ b/src/grid_impl/single_element_grid/geometry.rs @@ -1,49 +1,81 @@ //! Implementation of grid geometry -use crate::grid_impl::traits::Geometry; +use crate::grid_impl::common::{compute_jacobian, compute_normal_from_jacobian23, compute_point}; +use crate::grid_impl::traits::{Geometry, GeometryEvaluator}; +use crate::reference_cell; use bempp_element::element::CiarletElement; use bempp_traits::element::FiniteElement; use num::Float; +use rlst_common::types::Scalar; +use rlst_dense::{ + array::Array, + base_array::BaseArray, + data_container::VectorContainer, + rlst_array_from_slice2, rlst_dynamic_array4, + traits::{RandomAccessByRef, Shape, UnsafeRandomAccessByRef}, +}; /// Geometry of a serial grid -pub struct SerialSingleElementGeometry { +pub struct SerialSingleElementGeometry { dim: usize, index_map: Vec, - // TODO: change storage to rlst - coordinates: Vec, - cells: Vec, - element: CiarletElement, + pub(crate) coordinates: Array, 2>, 2>, + pub(crate) cells: Vec, + pub(crate) element: CiarletElement, midpoints: Vec>, diameters: Vec, volumes: Vec, + cell_indices: Vec, } -unsafe impl Sync for SerialSingleElementGeometry {} +unsafe impl Sync for SerialSingleElementGeometry {} -impl SerialSingleElementGeometry { +impl SerialSingleElementGeometry { pub fn new( - coordinates: Vec, - dim: usize, + coordinates: Array, 2>, 2>, cells_input: &[usize], - element: CiarletElement, + element: CiarletElement, ) -> Self { + let dim = coordinates.shape()[1]; let size = element.dim(); let ncells = cells_input.len() / size; let mut index_map = vec![0; ncells]; let mut cells = vec![]; - let mut midpoints = vec![]; + let mut midpoints = vec![vec![T::from(0.0).unwrap(); dim]; ncells]; let mut diameters = vec![]; let mut volumes = vec![]; - for (cell_i, i) in index_map.iter_mut().enumerate() { - *i = cell_i; - midpoints.push(vec![T::from(0.0).unwrap(); dim]); // TODO + let mut table = rlst_dynamic_array4!(T, element.tabulate_array_shape(0, 1)); + element.tabulate( + &rlst_array_from_slice2!( + T, + &reference_cell::midpoint(element.cell_type()), + [1, reference_cell::dim(element.cell_type())] + ), + 0, + &mut table, + ); + + let mut start = 0; + for (cell_i, index) in index_map.iter_mut().enumerate() { + *index = cell_i; + + for (i, v) in cells_input[start..start + size].iter().enumerate() { + let t = unsafe { *table.get_unchecked([0, 0, i, 0]) }; + for (j, component) in midpoints[cell_i].iter_mut().enumerate() { + *component += unsafe { *coordinates.get_unchecked([*v, j]) } * t; + } + } + diameters.push(T::from(0.0).unwrap()); // TODO volumes.push(T::from(0.0).unwrap()); // TODO + start += size; } cells.extend_from_slice(cells_input); + let cell_indices = (0..ncells).collect::>(); + Self { dim, index_map, @@ -53,14 +85,16 @@ impl SerialSingleElementGeometry { midpoints, diameters, volumes, + cell_indices, } } } -impl Geometry for SerialSingleElementGeometry { +impl Geometry for SerialSingleElementGeometry { type IndexType = usize; type T = T; - type Element = CiarletElement; + type Element = CiarletElement; + type Evaluator<'a> = GeometryEvaluatorSingleElement<'a, T> where Self: 'a; fn dim(&self) -> usize { self.dim @@ -71,15 +105,11 @@ impl Geometry for SerialSingleElementGeometry { } fn coordinate(&self, point_index: usize, coord_index: usize) -> Option<&Self::T> { - if coord_index < self.dim && point_index * self.dim < self.coordinates.len() { - Some(&self.coordinates[point_index * self.dim + coord_index]) - } else { - None - } + self.coordinates.get([point_index, coord_index]) } fn point_count(&self) -> usize { - self.coordinates.len() / self.dim + self.coordinates.shape()[0] } fn cell_points(&self, index: usize) -> Option<&[usize]> { @@ -107,15 +137,15 @@ impl Geometry for SerialSingleElementGeometry { 1 } fn element(&self, i: usize) -> Option<&Self::Element> { - if i < self.cell_count() { + if i == 0 { Some(&self.element) } else { None } } - fn cells(&self, i: usize) -> Option<&[usize]> { + fn cell_indices(&self, i: usize) -> Option<&[usize]> { if i == 0 { - Some(&self.cells) + Some(&self.cell_indices) } else { None } @@ -131,40 +161,154 @@ impl Geometry for SerialSingleElementGeometry { fn volume(&self, index: usize) -> Self::T { self.volumes[index] } + + fn get_evaluator<'a>(&'a self, points: &'a [Self::T]) -> GeometryEvaluatorSingleElement<'a, T> { + GeometryEvaluatorSingleElement::::new(self, points) + } +} + +pub struct GeometryEvaluatorSingleElement<'a, T: Float + Scalar> { + geometry: &'a SerialSingleElementGeometry, + tdim: usize, + table: Array, 4>, 4>, +} + +impl<'a, T: Float + Scalar> GeometryEvaluatorSingleElement<'a, T> { + fn new(geometry: &'a SerialSingleElementGeometry, points: &'a [T]) -> Self { + let tdim = reference_cell::dim(geometry.element.cell_type()); + assert_eq!(points.len() % tdim, 0); + let npoints = points.len() / tdim; + let rlst_points = rlst_array_from_slice2!(T, points, [tdim, npoints]); + + let mut table = rlst_dynamic_array4!(T, geometry.element.tabulate_array_shape(1, npoints)); + geometry.element.tabulate(&rlst_points, 1, &mut table); + Self { + geometry, + tdim, + table, + } + } +} + +impl<'a, T: Float + Scalar> GeometryEvaluator for GeometryEvaluatorSingleElement<'a, T> { + type T = T; + + fn point_count(&self) -> usize { + self.table.shape()[1] + } + + fn compute_point(&self, cell_index: usize, point_index: usize, point: &mut [T]) { + compute_point( + self.geometry, + self.table.view(), + cell_index, + point_index, + point, + ); + } + + fn compute_jacobian(&self, cell_index: usize, point_index: usize, jacobian: &mut [T]) { + compute_jacobian( + self.geometry, + self.table.view(), + self.tdim, + cell_index, + point_index, + jacobian, + ); + } + + fn compute_normal(&self, cell_index: usize, point_index: usize, normal: &mut [T]) { + let gdim = self.geometry.dim(); + let tdim = self.tdim; + assert_eq!(tdim, 2); + assert_eq!(tdim, gdim - 1); + + let mut jacobian = vec![T::from(0.0).unwrap(); gdim * tdim]; + self.compute_jacobian(cell_index, point_index, &mut jacobian[..]); + compute_normal_from_jacobian23(&jacobian, normal); + } } #[cfg(test)] mod test { use crate::grid_impl::single_element_grid::geometry::*; + use crate::types::ReferenceCellType; use approx::*; use bempp_element::element::{create_element, ElementFamily}; use bempp_traits::element::Continuity; + use rlst_dense::{ + rlst_dynamic_array2, + traits::{RandomAccessMut, RawAccess, RawAccessMut}, + }; - fn example_geometry() -> SerialSingleElementGeometry { + fn example_geometry_2d() -> SerialSingleElementGeometry { let p1triangle = create_element( ElementFamily::Lagrange, - bempp_element::cell::ReferenceCellType::Triangle, + ReferenceCellType::Triangle, 1, Continuity::Continuous, ); - SerialSingleElementGeometry::new( - vec![0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0], + let mut points = rlst_dynamic_array2!(f64, [4, 2]); + *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, 1]).unwrap() = 0.0; + *points.get_mut([2, 0]).unwrap() = 1.0; + *points.get_mut([2, 1]).unwrap() = 1.0; + *points.get_mut([3, 0]).unwrap() = 0.0; + *points.get_mut([3, 1]).unwrap() = 1.0; + SerialSingleElementGeometry::new(points, &[0, 1, 2, 0, 2, 3], p1triangle) + } + + fn example_geometry_3d() -> SerialSingleElementGeometry { + let p2triangle = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, 2, - &[0, 1, 2, 0, 2, 3], - p1triangle, - ) + Continuity::Continuous, + ); + let mut points = rlst_dynamic_array2!(f64, [9, 3]); + *points.get_mut([0, 0]).unwrap() = 0.0; + *points.get_mut([0, 1]).unwrap() = 0.0; + *points.get_mut([0, 2]).unwrap() = 0.0; + *points.get_mut([1, 0]).unwrap() = 0.5; + *points.get_mut([1, 1]).unwrap() = 0.0; + *points.get_mut([1, 2]).unwrap() = 0.5; + *points.get_mut([2, 0]).unwrap() = 1.0; + *points.get_mut([2, 1]).unwrap() = 0.0; + *points.get_mut([2, 2]).unwrap() = 0.0; + *points.get_mut([3, 0]).unwrap() = 0.0; + *points.get_mut([3, 1]).unwrap() = 0.5; + *points.get_mut([3, 2]).unwrap() = 0.0; + *points.get_mut([4, 0]).unwrap() = 0.5; + *points.get_mut([4, 1]).unwrap() = 0.5; + *points.get_mut([4, 2]).unwrap() = 0.0; + *points.get_mut([5, 0]).unwrap() = 1.0; + *points.get_mut([5, 1]).unwrap() = 0.5; + *points.get_mut([5, 2]).unwrap() = 0.0; + *points.get_mut([6, 0]).unwrap() = 0.0; + *points.get_mut([6, 1]).unwrap() = 1.0; + *points.get_mut([6, 2]).unwrap() = 0.0; + *points.get_mut([7, 0]).unwrap() = 0.5; + *points.get_mut([7, 1]).unwrap() = 1.0; + *points.get_mut([7, 2]).unwrap() = 0.0; + *points.get_mut([8, 0]).unwrap() = 1.0; + *points.get_mut([8, 1]).unwrap() = 1.0; + *points.get_mut([8, 2]).unwrap() = 0.0; + SerialSingleElementGeometry::new(points, &[0, 2, 8, 5, 4, 1, 0, 8, 6, 7, 3, 4], p2triangle) } #[test] fn test_counts() { - let g = example_geometry(); + let g = example_geometry_2d(); assert_eq!(g.point_count(), 4); assert_eq!(g.cell_count(), 2); } #[test] fn test_cell_points() { - let g = example_geometry(); + let g = example_geometry_2d(); for (cell_i, points) in [ vec![vec![0.0, 0.0], vec![1.0, 0.0], vec![1.0, 1.0]], vec![vec![0.0, 0.0], vec![1.0, 1.0], vec![0.0, 1.0]], @@ -175,9 +319,178 @@ mod test { let vs = g.cell_points(cell_i).unwrap(); for (p_i, point) in points.iter().enumerate() { for (c_i, coord) in point.iter().enumerate() { - assert_relative_eq!(*coord, *g.coordinate(vs[p_i], c_i).unwrap()); + assert_relative_eq!( + *coord, + *g.coordinate(vs[p_i], c_i).unwrap(), + epsilon = 1e-12 + ); } } } } + + fn triangle_points() -> Array, 2>, 2> { + let mut points = rlst_dynamic_array2!(f64, [2, 2]); + *points.get_mut([0, 0]).unwrap() = 0.2; + *points.get_mut([0, 1]).unwrap() = 0.5; + *points.get_mut([1, 0]).unwrap() = 0.6; + *points.get_mut([1, 1]).unwrap() = 0.1; + points + } + + #[test] + fn test_compute_point_2d() { + let g = example_geometry_2d(); + let points = triangle_points(); + + let evaluator = g.get_evaluator(points.data()); + let mut mapped_point = vec![0.0; 2]; + for (cell_i, points) in [ + vec![vec![0.7, 0.5], vec![0.7, 0.1]], + vec![vec![0.2, 0.7], vec![0.6, 0.7]], + ] + .iter() + .enumerate() + { + for (point_i, point) in points.iter().enumerate() { + evaluator.compute_point(cell_i, point_i, &mut mapped_point); + for (i, j) in mapped_point.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_compute_point_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut mapped_point = vec![0.0; 3]; + for (cell_i, points) in [ + vec![vec![0.7, 0.5, 0.12], vec![0.7, 0.1, 0.36]], + vec![vec![0.2, 0.7, 0.0], vec![0.6, 0.7, 0.0]], + ] + .iter() + .enumerate() + { + for (point_i, point) in points.iter().enumerate() { + evaluator.compute_point(cell_i, point_i, &mut mapped_point); + for (i, j) in mapped_point.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_compute_jacobian_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut computed_jacobian = rlst_dynamic_array2!(f64, [3, 2]); + for (cell_i, jacobians) in [ + vec![ + vec![vec![1.0, 1.0], vec![0.0, 1.0], vec![0.2, -0.4]], + vec![vec![1.0, 1.0], vec![0.0, 1.0], vec![-0.6, -1.2]], + ], + vec![ + vec![vec![1.0, 0.0], vec![1.0, 1.0], vec![0.0, 0.0]], + vec![vec![1.0, 0.0], vec![1.0, 1.0], vec![0.0, 0.0]], + ], + ] + .iter() + .enumerate() + { + for (point_i, jacobian) in jacobians.iter().enumerate() { + evaluator.compute_jacobian(cell_i, point_i, computed_jacobian.data_mut()); + for (i, row) in jacobian.iter().enumerate() { + for (j, entry) in row.iter().enumerate() { + assert_relative_eq!( + *entry, + *computed_jacobian.get([i, j]).unwrap(), + epsilon = 1e-12 + ); + } + } + } + } + } + #[test] + fn test_compute_normal_3d() { + let g = example_geometry_3d(); + let points = triangle_points(); + let evaluator = g.get_evaluator(points.data()); + + let mut computed_normal = vec![0.0; 3]; + for (cell_i, normals) in [ + vec![ + vec![ + -0.2 / f64::sqrt(1.4), + 0.6 / f64::sqrt(1.4), + 1.0 / f64::sqrt(1.4), + ], + vec![ + 0.6 / f64::sqrt(1.72), + 0.6 / f64::sqrt(1.72), + 1.0 / f64::sqrt(1.72), + ], + ], + vec![vec![0.0, 0.0, 1.0], vec![0.0, 0.0, 1.0]], + ] + .iter() + .enumerate() + { + for (point_i, normal) in normals.iter().enumerate() { + evaluator.compute_normal(cell_i, point_i, &mut computed_normal); + assert_relative_eq!( + computed_normal[0] * computed_normal[0] + + computed_normal[1] * computed_normal[1] + + computed_normal[2] * computed_normal[2], + 1.0, + epsilon = 1e-12 + ); + for (i, j) in computed_normal.iter().zip(normal) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + } + + #[test] + fn test_midpoint_2d() { + let g = example_geometry_2d(); + + let mut midpoint = vec![0.0; 2]; + for (cell_i, point) in [vec![2.0 / 3.0, 1.0 / 3.0], vec![1.0 / 3.0, 2.0 / 3.0]] + .iter() + .enumerate() + { + g.midpoint(cell_i, &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } + + #[test] + fn test_midpoint_3d() { + let g = example_geometry_3d(); + + let mut midpoint = vec![0.0; 3]; + for (cell_i, point) in [ + vec![2.0 / 3.0, 1.0 / 3.0, 2.0 / 9.0], + vec![1.0 / 3.0, 2.0 / 3.0, 0.0], + ] + .iter() + .enumerate() + { + g.midpoint(cell_i, &mut midpoint); + for (i, j) in midpoint.iter().zip(point) { + assert_relative_eq!(*i, *j, epsilon = 1e-12); + } + } + } } diff --git a/src/grid_impl/single_element_grid/grid.rs b/src/grid_impl/single_element_grid/grid.rs index c96af4f..5da00d3 100644 --- a/src/grid_impl/single_element_grid/grid.rs +++ b/src/grid_impl/single_element_grid/grid.rs @@ -6,37 +6,32 @@ use crate::grid_impl::single_element_grid::{ use crate::grid_impl::traits::Grid; use crate::reference_cell; use crate::reference_cell::ReferenceCellType; -use bempp_element::element::{create_element, ElementFamily}; +use bempp_element::element::{create_element, ElementFamily, Inverse}; use bempp_traits::element::{Continuity, FiniteElement}; +use log::warn; use num::Float; +use rlst_common::types::Scalar; +use rlst_dense::{array::Array, base_array::BaseArray, data_container::VectorContainer}; /// A serial grid -pub struct SerialSingleElementGrid { +pub struct SerialSingleElementGrid { topology: SerialSingleElementTopology, geometry: SerialSingleElementGeometry, } -impl SerialSingleElementGrid { +impl SerialSingleElementGrid { pub fn new( - points: Vec, - gdim: usize, + points: Array, 2>, 2>, cells: &[usize], cell_type: ReferenceCellType, cell_degree: usize, ) -> Self { - let element = create_element( + if cell_type == ReferenceCellType::Triangle && cell_degree == 1 { + warn!("Creating a single element grid with a P1 triangle. Using a FlatTriangleGrid would be faster."); + } + let element = create_element::( ElementFamily::Lagrange, - // TODO: remove this match once bempp-rs and grid-rs use the same ReferenceCellType - match cell_type { - ReferenceCellType::Interval => bempp_element::cell::ReferenceCellType::Interval, - ReferenceCellType::Triangle => bempp_element::cell::ReferenceCellType::Triangle, - ReferenceCellType::Quadrilateral => { - bempp_element::cell::ReferenceCellType::Quadrilateral - } - _ => { - panic!("Unsupported cell type: {:?}", cell_type); - } - }, + cell_type, cell_degree, Continuity::Continuous, ); @@ -55,14 +50,14 @@ impl SerialSingleElementGrid { let topology = SerialSingleElementTopology::new(&cell_vertices, cell_type); // Create the geometry - let geometry = SerialSingleElementGeometry::::new(points, gdim, cells, element); + let geometry = SerialSingleElementGeometry::::new(points, cells, element); Self { topology, geometry } } } /// A grid -impl Grid for SerialSingleElementGrid { +impl Grid for SerialSingleElementGrid { type T = T; /// The type that implements [Topology] diff --git a/src/grid_impl/single_element_grid/topology.rs b/src/grid_impl/single_element_grid/topology.rs index 5b29e5b..1f9bef1 100644 --- a/src/grid_impl/single_element_grid/topology.rs +++ b/src/grid_impl/single_element_grid/topology.rs @@ -3,6 +3,7 @@ use crate::grid_impl::traits::{Ownership, Topology}; use crate::reference_cell; use crate::reference_cell::ReferenceCellType; +use crate::types::CellLocalIndexPair; fn all_equal(a: &[T], b: &[T]) -> bool { if a.len() != b.len() { @@ -25,9 +26,9 @@ fn all_in(a: &[T], b: &[T]) -> bool { pub struct SerialSingleElementTopology { dim: usize, index_map: Vec, - cells: Vec>, // TODO: use 2D array - connectivity: Vec>>>, - cell_connectivity: Vec>>, + entities_to_vertices: Vec>>, + cells_to_entities: Vec>>, + entities_to_cells: Vec>>>, entity_types: Vec, } @@ -47,13 +48,13 @@ impl SerialSingleElementTopology { .filter(|t| !t.is_empty()) .map(|t| t[0]) .collect::>(); - let mut cells = vec![]; - let mut connectivity = vec![vec![vec![]; dim + 1]; dim + 1]; - let mut cell_connectivity = vec![vec![]; dim + 1]; - // dim0 = dim, dim1 = 0 - let mut cty = vec![]; - let mut cell_cty = vec![vec![]; dim + 1]; + let mut cells_to_entities = vec![vec![vec![]; ncells]; dim + 1]; + let mut entities_to_cells = vec![vec![]; dim + 1]; + let mut entities_to_vertices = vec![vec![]; dim]; + + entities_to_cells[dim] = vec![vec![]; ncells]; + let mut start = 0; for (cell_i, i) in index_map.iter_mut().enumerate() { let cell = &cells_input[start..start + size]; @@ -61,130 +62,61 @@ impl SerialSingleElementTopology { let mut row = vec![]; for v in cell { if !vertices.contains(v) { + entities_to_cells[0].push(vec![]); vertices.push(*v); } row.push(vertices.iter().position(|&r| r == *v).unwrap()); } - cell_cty[0].extend_from_slice(&row); - cells.push(row.clone()); - cty.push(row); + + for (local_index, v) in row.iter().enumerate() { + entities_to_cells[0][*v].push(CellLocalIndexPair::new(cell_i, local_index)); + } + entities_to_cells[dim][cell_i] = vec![CellLocalIndexPair::new(cell_i, 0)]; + + cells_to_entities[0][cell_i] = row; + cells_to_entities[dim][cell_i] = vec![cell_i]; start += size; } - connectivity[dim][0] = cty; - cell_connectivity[dim] = cell_cty; - // dim1 == 0 - for dim0 in 0..dim { - let mut cty: Vec> = vec![]; - let ref_conn = &reference_cell::connectivity(cell_type)[dim0]; - for cell in &connectivity[dim][0] { - for rc in ref_conn { + + for i in 0..vertices.len() { + entities_to_vertices[0].push(vec![i]); + } + for d in 1..dim { + let mut c_to_e = vec![]; + let ref_conn = &reference_cell::connectivity(cell_type)[d]; + for (cell_i, cell) in cells_to_entities[0].iter().enumerate() { + let mut entity_ids = vec![]; + for (local_index, rc) in ref_conn.iter().enumerate() { let vertices = rc[0].iter().map(|x| cell[*x]).collect::>(); let mut found = false; - for entity in &cty { + for (entity_index, entity) in entities_to_vertices[d].iter().enumerate() { if all_equal(entity, &vertices) { + entity_ids.push(entity_index); + entities_to_cells[d][entity_index] + .push(CellLocalIndexPair::new(cell_i, local_index)); found = true; break; } } if !found { - cty.push(vertices); + entity_ids.push(entities_to_vertices[d].len()); + entities_to_cells[d] + .push(vec![CellLocalIndexPair::new(cell_i, local_index)]); + entities_to_vertices[d].push(vertices); } } + c_to_e.push(entity_ids); } - connectivity[dim0][0] = cty; - } - // dim0 == dim1 == 0 - let mut nvertices = 0; - let mut cty = vec![]; - for cell in &connectivity[dim][0] { - nvertices = std::cmp::max(nvertices, 1 + cell.iter().max().unwrap()); - } - for i in 0..nvertices { - cty.push(vec![i]); - } - connectivity[0][0] = cty; - // dim0 == dim1 - for (dim0, c) in connectivity.iter_mut().enumerate() { - let mut cty = vec![]; - for i in 0..c[0].len() { - cty.push(vec![i]); - } - c[dim0] = cty; - } - // dim0 == dim1 == dim - let mut cty = vec![]; - for i in 0..connectivity[dim][0].len() { - cty.push(i); + cells_to_entities[d] = c_to_e; } - cell_connectivity[dim][dim] = cty; - // dim0 == dim - for dim1 in 1..dim { - let mut cty = vec![]; - let mut cell_cty = vec![]; - let entities0 = &connectivity[dim][0]; - let ref_conn = &reference_cell::connectivity(cell_type)[dim1]; - let entities1 = &connectivity[dim1][0]; - - for entity0 in entities0 { - let mut row = vec![]; - for rc in ref_conn { - let vertices = rc[0].iter().map(|x| entity0[*x]).collect::>(); - for (j, entity1) in entities1.iter().enumerate() { - if all_equal(&vertices, entity1) { - row.push(j); - break; - } - } - } - cell_cty.extend_from_slice(&row); - cty.push(row); - } - cell_connectivity[dim][dim1] = cell_cty; - connectivity[dim][dim1] = cty; - } - // dim1 < dim0 - for (dim0, etype0) in entity_types.iter().enumerate().take(dim + 1).skip(2) { - for dim1 in 1..dim0 { - let mut cty = vec![]; - let entities0 = &connectivity[dim0][0]; - let ref_conn = &reference_cell::connectivity(*etype0)[dim1]; - let entities1 = &connectivity[dim1][0]; - for entity0 in entities0 { - let mut row = vec![]; - for rc in ref_conn { - let vertices = rc[0].iter().map(|x| entity0[*x]).collect::>(); - for (j, entity1) in entities1.iter().enumerate() { - if all_equal(&vertices, entity1) { - row.push(j); - break; - } - } - } - cty.push(row); - } - connectivity[dim0][dim1] = cty; - } - } - // dim1 > dim0 - for dim0 in 0..dim { - for dim1 in 1..dim + 1 { - let mut data = vec![vec![]; connectivity[dim0][0].len()]; - for (i, row) in connectivity[dim1][dim0].iter().enumerate() { - for v in row { - data[*v].push(i); - } - } - connectivity[dim0][dim1] = data; - } - } Self { dim, index_map, - cells, - connectivity, - cell_connectivity, + entities_to_vertices, + cells_to_entities, + entities_to_cells, entity_types, } } @@ -201,7 +133,7 @@ impl Topology for SerialSingleElementTopology { } fn entity_count(&self, etype: ReferenceCellType) -> usize { if self.entity_types.contains(&etype) { - self.connectivity[reference_cell::dim(etype)][0].len() + self.entities_to_cells[reference_cell::dim(etype)].len() } else { 0 } @@ -210,14 +142,14 @@ impl Topology for SerialSingleElementTopology { self.entity_count(self.entity_types[dim]) } fn cell(&self, index: Self::IndexType) -> Option<&[usize]> { - if index < self.cells.len() { - Some(&self.cells[index]) + if index < self.cells_to_entities[self.dim].len() { + Some(&self.cells_to_entities[self.dim][index]) } else { None } } fn cell_type(&self, index: Self::IndexType) -> Option { - if index < self.cells.len() { + if index < self.cells_to_entities[self.dim].len() { Some(self.entity_types[self.dim]) } else { None @@ -228,33 +160,36 @@ impl Topology for SerialSingleElementTopology { &self.entity_types[dim..dim + 1] } - fn cell_entities(&self, cell_type: ReferenceCellType, dim: usize) -> Option<&[usize]> { - if self.entity_types.contains(&cell_type) - && dim < self.cell_connectivity[reference_cell::dim(cell_type)].len() - { - Some(&self.cell_connectivity[reference_cell::dim(cell_type)][dim]) + fn entity_ownership(&self, _dim: usize, _index: Self::IndexType) -> Ownership { + Ownership::Owned + } + fn cell_to_entities(&self, index: Self::IndexType, dim: usize) -> Option<&[Self::IndexType]> { + if dim <= self.dim && index < self.cells_to_entities[dim].len() { + Some(&self.cells_to_entities[dim][index]) } else { None } } - - fn connectivity(&self, dim0: usize, index: usize, dim1: usize) -> Option<&[Self::IndexType]> { - if dim0 < self.connectivity.len() - && dim1 < self.connectivity[dim0].len() - && index < self.connectivity[dim0][dim1].len() - { - Some(&self.connectivity[dim0][dim1][index]) + fn entity_to_cells( + &self, + dim: usize, + index: Self::IndexType, + ) -> Option<&[CellLocalIndexPair]> { + if dim <= self.dim && index < self.entities_to_cells[dim].len() { + Some(&self.entities_to_cells[dim][index]) } else { None } } - fn entity_ownership(&self, _dim: usize, _index: Self::IndexType) -> Ownership { - Ownership::Owned - } - - fn extract_index(&self, index: usize) -> usize { - index + fn entity_vertices(&self, dim: usize, index: Self::IndexType) -> Option<&[Self::IndexType]> { + if dim == self.dim { + self.cell_to_entities(index, 0) + } else if dim < self.dim && index < self.entities_to_vertices[dim].len() { + Some(&self.entities_to_vertices[dim][index]) + } else { + None + } } } @@ -278,136 +213,85 @@ mod test { #[test] fn test_cell_entities_points() { let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 0).unwrap(); - assert_eq!(c.len(), 6); - // cell 0 - assert_eq!(c[0], 0); - assert_eq!(c[1], 1); - assert_eq!(c[2], 2); - // cell 1 - assert_eq!(c[3], 2); - assert_eq!(c[4], 1); - assert_eq!(c[5], 3); + for (i, vertices) in [[0, 1, 2], [2, 1, 3]].iter().enumerate() { + let c = t.cell_to_entities(i, 0).unwrap(); + assert_eq!(c.len(), 3); + assert_eq!(c, vertices); + } } #[test] fn test_cell_entities_intervals() { let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 1).unwrap(); - assert_eq!(c.len(), 6); - // cell 0 - assert_eq!(c[0], 0); - assert_eq!(c[1], 1); - assert_eq!(c[2], 2); - // cell 1 - assert_eq!(c[3], 3); - assert_eq!(c[4], 4); - assert_eq!(c[5], 0); + for (i, edges) in [[0, 1, 2], [3, 4, 0]].iter().enumerate() { + let c = t.cell_to_entities(i, 1).unwrap(); + assert_eq!(c.len(), 3); + assert_eq!(c, edges); + } } + #[test] fn test_cell_entities_triangles() { let t = example_topology(); - let c = t.cell_entities(ReferenceCellType::Triangle, 2).unwrap(); - assert_eq!(c.len(), 2); - // cell 0 - assert_eq!(c[0], 0); - // cell 1 - assert_eq!(c[1], 1); - } - #[test] - fn test_vertex_connectivity() { - let t = example_topology(); - - for (id, vertices) in [vec![0], vec![1], vec![2], vec![3]].iter().enumerate() { - let c = t.connectivity(0, id, 0).unwrap(); - for (i, j) in c.iter().zip(vertices) { - assert_eq!(*i, *j); - } - } - - for (id, edges) in [vec![1, 2], vec![0, 2, 3], vec![0, 1, 4], vec![3, 4]] - .iter() - .enumerate() - { - let c = t.connectivity(0, id, 1).unwrap(); - for (i, j) in c.iter().zip(edges) { - assert_eq!(*i, *j); - } - } - - for (id, faces) in [vec![0], vec![0, 1], vec![0, 1], vec![1]] - .iter() - .enumerate() - { - let c = t.connectivity(0, id, 2).unwrap(); - for (i, j) in c.iter().zip(faces) { - assert_eq!(*i, *j); - } + for i in 0..2 { + let c = t.cell_to_entities(i, 2).unwrap(); + assert_eq!(c.len(), 1); + assert_eq!(c[0], i); } } + #[test] - fn test_edge_connectivity() { + fn test_entities_to_cells_points() { let t = example_topology(); + let c_to_e = (0..t.entity_count(ReferenceCellType::Triangle)) + .map(|i| t.cell_to_entities(i, 0).unwrap()) + .collect::>(); + let e_to_c = (0..t.entity_count(ReferenceCellType::Point)) + .map(|i| { + t.entity_to_cells(0, i) + .unwrap() + .iter() + .map(|x| x.cell) + .collect::>() + }) + .collect::>(); - for (id, vertices) in [vec![1, 2], vec![0, 2], vec![0, 1], vec![1, 3], vec![2, 3]] - .iter() - .enumerate() - { - let c = t.connectivity(1, id, 0).unwrap(); - for (i, j) in c.iter().zip(vertices) { - assert_eq!(*i, *j); + for (i, cell) in c_to_e.iter().enumerate() { + for v in *cell { + assert!(e_to_c[*v].contains(&i)); } } - - for (id, edges) in [vec![0], vec![1], vec![2], vec![3], vec![4]] - .iter() - .enumerate() - { - let c = t.connectivity(1, id, 1).unwrap(); - for (i, j) in c.iter().zip(edges) { - assert_eq!(*i, *j); - } - } - - for (id, faces) in [vec![0, 1], vec![0], vec![0], vec![1], vec![1]] - .iter() - .enumerate() - { - let c = t.connectivity(1, id, 2).unwrap(); - for (i, j) in c.iter().zip(faces) { - assert_eq!(*i, *j); + for (i, cells) in e_to_c.iter().enumerate() { + for c in cells { + assert!(c_to_e[*c].contains(&i)); } } } #[test] - fn test_cell_entities_vs_connectivity() { + fn test_entities_to_cells_edges() { let t = example_topology(); - for cell_type in t.entity_types(t.dim()) { - for dim in 0..t.dim() + 1 { - let ce = t.cell_entities(*cell_type, dim).unwrap(); - let n = reference_cell::entity_counts(*cell_type)[dim]; - for i in 0..ce.len() / n { - // TODO: entity_count instead? - println!("{cell_type:?} {dim} {i}"); - let con = t.connectivity(t.dim(), i, dim).unwrap(); - println!( - "{cell_type:?} {dim} {i} -> {:?} {:?}", - con, - &ce[n * i..n * (i + 1)] - ); - } + let c_to_e = (0..t.entity_count(ReferenceCellType::Triangle)) + .map(|i| t.cell_to_entities(i, 1).unwrap()) + .collect::>(); + let e_to_c = (0..t.entity_count(ReferenceCellType::Interval)) + .map(|i| { + t.entity_to_cells(1, i) + .unwrap() + .iter() + .map(|x| x.cell) + .collect::>() + }) + .collect::>(); + + for (i, cell) in c_to_e.iter().enumerate() { + for v in *cell { + assert!(e_to_c[*v].contains(&i)); } } - for cell_type in t.entity_types(t.dim()) { - for dim in 0..t.dim() + 1 { - let ce = t.cell_entities(*cell_type, dim).unwrap(); - let n = reference_cell::entity_counts(*cell_type)[dim]; - for i in 0..ce.len() / n { - // TODO: entity_count instead? - let con = t.connectivity(t.dim(), i, dim).unwrap(); - assert_eq!(con, &ce[n * i..n * (i + 1)]); - } + for (i, cells) in e_to_c.iter().enumerate() { + for c in cells { + assert!(c_to_e[*c].contains(&i)); } } } diff --git a/src/grid_impl/traits.rs b/src/grid_impl/traits.rs index 91dfa1a..44ec35f 100644 --- a/src/grid_impl/traits.rs +++ b/src/grid_impl/traits.rs @@ -1,6 +1,7 @@ //! Traits used in the implementation of a grid use crate::reference_cell::ReferenceCellType; +use crate::types::CellLocalIndexPair; use bempp_traits::element::FiniteElement; use num::Float; @@ -37,23 +38,21 @@ pub trait Topology { /// All entity types of the given dimension that are included in the grid fn entity_types(&self, dim: usize) -> &[ReferenceCellType]; - /// Get the indices of entities and types of entities that are connected to each cell of type `cell_type` - fn cell_entities(&self, cell_type: ReferenceCellType, dim: usize) - -> Option<&[Self::IndexType]>; + /// Get the indices of entities of dimension `dim` that are connected to the cell with index `index` + fn cell_to_entities(&self, index: Self::IndexType, dim: usize) -> Option<&[Self::IndexType]>; - /// Get the indices and types of entities of dimension `dim1` that are connected to the entity of dimension `dim0` with index `index` - fn connectivity( + /// Get the indices of entities of cell that are connected to the entity with dimension `dim` and index `index` + fn entity_to_cells( &self, - dim0: usize, + dim: usize, index: Self::IndexType, - dim1: usize, - ) -> Option<&[Self::IndexType]>; + ) -> Option<&[CellLocalIndexPair]>; + + /// Get the indices of the vertices that are connect to theentity with dimension `dim` and index `index` + fn entity_vertices(&self, dim: usize, index: Self::IndexType) -> Option<&[Self::IndexType]>; /// Get the ownership of a mesh entity fn entity_ownership(&self, dim: usize, index: Self::IndexType) -> Ownership; - - /// Extract a flat index from an IndexType - fn extract_index(&self, index: Self::IndexType) -> usize; } /// The geometry of a grid @@ -63,6 +62,9 @@ pub trait Geometry { type IndexType: std::fmt::Debug + Eq + Copy; type T: Float; type Element: FiniteElement; + type Evaluator<'a>: GeometryEvaluator + where + Self: 'a; /// The geometric dimension fn dim(&self) -> usize; @@ -90,7 +92,7 @@ pub trait Geometry { /// Get the `i`th element fn element(&self, i: usize) -> Option<&Self::Element>; /// Get the cells associated with the `i`th element - fn cells(&self, i: usize) -> Option<&[usize]>; + fn cell_indices(&self, i: usize) -> Option<&[Self::IndexType]>; // ... or would it be better to replace the above 3 functions with an Iter? // fn elements_and_cells(&self) -> std::slice::Iter<'_, (&Self::Element, &[usize])>; @@ -103,6 +105,28 @@ pub trait Geometry { /// Volume of a cell fn volume(&self, index: Self::IndexType) -> Self::T; + + /// Get the geometry evaluator for the given points + fn get_evaluator<'a>(&'a self, points: &'a [Self::T]) -> Self::Evaluator<'a>; +} + +/// Geometry evaluator +/// +/// A geometry evaluator can compute points and jacobians on physical cells +pub trait GeometryEvaluator { + type T: Float; + + /// The number of points on the reference cell used by this evaluator + fn point_count(&self) -> usize; + + /// Compute a point on a physical cell + fn compute_point(&self, cell_index: usize, point_index: usize, point: &mut [Self::T]); + + /// Compute a jacobian on a physical cell + fn compute_jacobian(&self, cell_index: usize, point_index: usize, jacobian: &mut [Self::T]); + + /// Compute a normal on a physical cell + fn compute_normal(&self, cell_index: usize, point_index: usize, normal: &mut [Self::T]); } /// A grid diff --git a/src/reference_cell.rs b/src/reference_cell.rs index b50f377..10cf58e 100644 --- a/src/reference_cell.rs +++ b/src/reference_cell.rs @@ -56,6 +56,26 @@ pub fn vertices(cell: ReferenceCellType) -> Vec { } } +/// The midpoint of the cell +pub fn midpoint(cell: ReferenceCellType) -> Vec { + let half = T::from(0.5).unwrap(); + let third = T::from(1.0).unwrap() / T::from(3.0).unwrap(); + match cell { + ReferenceCellType::Point => vec![], + ReferenceCellType::Interval => vec![half], + ReferenceCellType::Triangle => vec![third; 2], + ReferenceCellType::Quadrilateral => vec![half; 2], + ReferenceCellType::Tetrahedron => vec![T::from(1.0).unwrap() / T::from(6.0).unwrap(); 3], + ReferenceCellType::Hexahedron => vec![half; 3], + ReferenceCellType::Prism => vec![third, third, half], + ReferenceCellType::Pyramid => vec![ + T::from(0.4).unwrap(), + T::from(0.4).unwrap(), + T::from(0.2).unwrap(), + ], + } +} + /// The edges of the reference cell pub fn edges(cell: ReferenceCellType) -> Vec> { match cell { diff --git a/src/traits/grid.rs b/src/traits/grid.rs index e1f2a9b..843ed7d 100644 --- a/src/traits/grid.rs +++ b/src/traits/grid.rs @@ -3,13 +3,14 @@ use crate::traits::cell::CellType; use crate::traits::point::PointType; use crate::types::cell_iterator::CellIterator; -use crate::types::vertex_iterator::PointIterator; +use crate::types::point_iterator::PointIterator; use crate::types::{CellLocalIndexPair, Float}; use super::ReferenceMapType; pub trait GridType: std::marker::Sized { type T: Float; + type IndexType: std::fmt::Debug + Eq + Copy; type Point<'a>: PointType where @@ -23,13 +24,6 @@ pub trait GridType: std::marker::Sized { where Self: 'a; - type ReferenceMapIterator<'a, Iter: std::iter::Iterator>: std::iter::Iterator< - Item = Self::ReferenceMap<'a>, - > - where - Self: 'a, - Iter: 'a; - fn number_of_vertices(&self) -> usize; fn number_of_points(&self) -> usize; @@ -71,20 +65,14 @@ pub trait GridType: std::marker::Sized { fn reference_to_physical_map<'a>( &'a self, reference_points: &'a [Self::T], - cell_index: usize, ) -> Self::ReferenceMap<'a>; - fn iter_reference_to_physical_map<'a, Iter: std::iter::Iterator + 'a>( - &'a self, - reference_points: &'a [Self::T], - iter: Iter, - ) -> Self::ReferenceMapIterator<'a, Iter> - where - Self: 'a; - - fn point_to_cells(&self, point_index: usize) -> &[CellLocalIndexPair]; + fn vertex_to_cells( + &self, + vertex_index: Self::IndexType, + ) -> &[CellLocalIndexPair]; - fn edge_to_cells(&self, edge_index: usize) -> &[CellLocalIndexPair]; + fn edge_to_cells(&self, edge_index: Self::IndexType) -> &[CellLocalIndexPair]; - fn face_to_cells(&self, face_index: usize) -> &[CellLocalIndexPair]; + fn face_to_cells(&self, face_index: Self::IndexType) -> &[CellLocalIndexPair]; } diff --git a/src/traits/reference_map.rs b/src/traits/reference_map.rs index 46b4313..e2f604e 100644 --- a/src/traits/reference_map.rs +++ b/src/traits/reference_map.rs @@ -15,13 +15,28 @@ pub trait ReferenceMapType { /// Defines an iterator that returns a slice with the value of the /// physical point for each reference point. - fn reference_to_physical(&self, point_index: usize, value: &mut [::T]); + fn reference_to_physical( + &self, + cell_index: usize, + point_index: usize, + value: &mut [::T], + ); /// Defines an iterator that returns a slice with the value of the /// Jacobian at the physical point for each reference point. - fn jacobian(&self, point_index: usize, value: &mut [::T]); + fn jacobian( + &self, + cell_index: usize, + point_index: usize, + value: &mut [::T], + ); /// Defines an iterator that returns a slice with the normal direction /// at each point. - fn normal(&self, point_index: usize, value: &mut [::T]); + fn normal( + &self, + cell_index: usize, + point_index: usize, + value: &mut [::T], + ); } diff --git a/src/types.rs b/src/types.rs index 686c3ca..397ee85 100644 --- a/src/types.rs +++ b/src/types.rs @@ -1,29 +1,20 @@ //! Type definitions pub mod cell_iterator; -pub mod vertex_iterator; +pub mod point_iterator; pub use num::Float; -#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] -pub enum ReferenceCellType { - Point, - Interval, - Triangle, - Quadrilateral, - Tetrahedron, - Hexahedron, - Prism, - Pyramid, -} +pub use bempp_element::cell::ReferenceCellType; -pub struct CellLocalIndexPair { - pub cell: usize, +#[derive(Debug, Clone)] +pub struct CellLocalIndexPair { + pub cell: IndexType, pub local_index: usize, } -impl CellLocalIndexPair { - pub fn new(cell: usize, local_index: usize) -> Self { +impl CellLocalIndexPair { + pub fn new(cell: IndexType, local_index: usize) -> Self { Self { cell, local_index } } } diff --git a/src/types/vertex_iterator.rs b/src/types/point_iterator.rs similarity index 100% rename from src/types/vertex_iterator.rs rename to src/types/point_iterator.rs