diff --git a/src/grid/single_element_grid/builder.rs b/src/grid/single_element_grid/builder.rs index 761ae40..6496216 100644 --- a/src/grid/single_element_grid/builder.rs +++ b/src/grid/single_element_grid/builder.rs @@ -22,8 +22,8 @@ pub struct SerialSingleElementGridBuilder, cells: Vec, point_indices_to_ids: Vec, - cell_indices_to_ids: Vec, point_ids_to_indices: HashMap, + cell_indices_to_ids: Vec, cell_ids_to_indices: HashMap, } @@ -51,8 +51,8 @@ where points: vec![], cells: vec![], point_indices_to_ids: vec![], - cell_indices_to_ids: vec![], point_ids_to_indices: HashMap::new(), + cell_indices_to_ids: vec![], cell_ids_to_indices: HashMap::new(), } } @@ -71,8 +71,8 @@ where points: Vec::with_capacity(npoints * Self::GDIM), cells: Vec::with_capacity(ncells * points_per_cell), point_indices_to_ids: Vec::with_capacity(npoints), - cell_indices_to_ids: Vec::with_capacity(ncells), point_ids_to_indices: HashMap::new(), + cell_indices_to_ids: Vec::with_capacity(ncells), cell_ids_to_indices: HashMap::new(), } } @@ -109,6 +109,10 @@ where &self.cells, self.element_data.0, self.element_data.1, + self.point_indices_to_ids, + self.point_ids_to_indices, + self.cell_indices_to_ids, + self.cell_ids_to_indices, ) } } diff --git a/src/grid/single_element_grid/geometry.rs b/src/grid/single_element_grid/geometry.rs index c99796d..3e0eb47 100644 --- a/src/grid/single_element_grid/geometry.rs +++ b/src/grid/single_element_grid/geometry.rs @@ -14,6 +14,7 @@ use rlst_dense::{ rlst_array_from_slice2, rlst_dynamic_array4, traits::{RandomAccessByRef, Shape, UnsafeRandomAccessByRef}, }; +use std::collections::HashMap; /// Geometry of a serial grid pub struct SerialSingleElementGeometry { @@ -26,6 +27,10 @@ pub struct SerialSingleElementGeometry { diameters: Vec, volumes: Vec, cell_indices: Vec, + point_indices_to_ids: Vec, + point_ids_to_indices: HashMap, + cell_indices_to_ids: Vec, + cell_ids_to_indices: HashMap, } unsafe impl Sync for SerialSingleElementGeometry {} @@ -35,6 +40,10 @@ impl SerialSingleElementGeometry { coordinates: Array, 2>, 2>, cells_input: &[usize], element: CiarletElement, + point_indices_to_ids: Vec, + point_ids_to_indices: HashMap, + cell_indices_to_ids: Vec, + cell_ids_to_indices: HashMap, ) -> Self { let dim = coordinates.shape()[1]; let size = element.dim(); @@ -86,6 +95,10 @@ impl SerialSingleElementGeometry { diameters, volumes, cell_indices, + point_indices_to_ids, + point_ids_to_indices, + cell_indices_to_ids, + cell_ids_to_indices, } } } @@ -167,18 +180,16 @@ impl Geometry for SerialSingleElementGeometry { } fn point_index_to_id(&self, index: usize) -> usize { - panic!(); + self.point_indices_to_ids[index] } - fn cell_index_to_id(&self, index: usize) -> usize { - panic!(); + self.cell_indices_to_ids[index] } - fn point_id_to_index(&self, id: usize) -> usize { - panic!(); + self.point_ids_to_indices[&id] } fn cell_id_to_index(&self, id: usize) -> usize { - panic!(); + self.cell_ids_to_indices[&id] } } @@ -273,7 +284,15 @@ mod test { *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) + SerialSingleElementGeometry::new( + points, + &[0, 1, 2, 0, 2, 3], + p1triangle, + vec![0, 1, 2, 3], + HashMap::from([(0, 0), (1, 1), (2, 2), (3, 3)]), + vec![0, 1], + HashMap::from([(0, 0), (1, 1)]), + ) } fn example_geometry_3d() -> SerialSingleElementGeometry { @@ -311,7 +330,25 @@ mod test { *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) + SerialSingleElementGeometry::new( + points, + &[0, 2, 8, 5, 4, 1, 0, 8, 6, 7, 3, 4], + p2triangle, + vec![0, 1, 2, 3, 4, 5, 6, 7, 8], + HashMap::from([ + (0, 0), + (1, 1), + (2, 2), + (3, 3), + (4, 4), + (5, 5), + (6, 6), + (7, 7), + (8, 8), + ]), + vec![0, 1], + HashMap::from([(0, 0), (1, 1)]), + ) } #[test] diff --git a/src/grid/single_element_grid/grid.rs b/src/grid/single_element_grid/grid.rs index 1e9a44a..720a011 100644 --- a/src/grid/single_element_grid/grid.rs +++ b/src/grid/single_element_grid/grid.rs @@ -17,6 +17,7 @@ use rlst_dense::{ data_container::VectorContainer, traits::MatrixInverse, }; +use std::collections::HashMap; /// A serial grid pub struct SerialSingleElementGrid { @@ -28,11 +29,17 @@ impl SerialSingleElementGrid where for<'a> Array, 2>, 2>, 2>: MatrixInverse, { + + #[allow(clippy::too_many_arguments)] pub fn new( points: Array, 2>, 2>, cells: &[usize], cell_type: ReferenceCellType, cell_degree: usize, + point_indices_to_ids: Vec, + point_ids_to_indices: HashMap, + cell_indices_to_ids: Vec, + cell_ids_to_indices: HashMap, ) -> Self { if cell_type == ReferenceCellType::Triangle && cell_degree == 1 { warn!("Creating a single element grid with a P1 triangle. Using a FlatTriangleGrid would be faster."); @@ -54,8 +61,21 @@ where start += npoints; } - let topology = SerialSingleElementTopology::new(&cell_vertices, cell_type); - let geometry = SerialSingleElementGeometry::::new(points, cells, element); + let topology = SerialSingleElementTopology::new( + &cell_vertices, + cell_type, + &point_indices_to_ids, + &cell_indices_to_ids, + ); + let geometry = SerialSingleElementGeometry::::new( + points, + cells, + element, + point_indices_to_ids, + point_ids_to_indices, + cell_indices_to_ids, + cell_ids_to_indices, + ); Self { topology, geometry } } } diff --git a/src/grid/single_element_grid/topology.rs b/src/grid/single_element_grid/topology.rs index 9b13853..c115a56 100644 --- a/src/grid/single_element_grid/topology.rs +++ b/src/grid/single_element_grid/topology.rs @@ -4,6 +4,7 @@ use crate::grid::traits::{Ownership, Topology}; use crate::reference_cell; use crate::reference_cell::ReferenceCellType; use crate::types::CellLocalIndexPair; +use std::collections::HashMap; fn all_equal(a: &[T], b: &[T]) -> bool { if a.len() != b.len() { @@ -30,15 +31,29 @@ pub struct SerialSingleElementTopology { cells_to_entities: Vec>>, entities_to_cells: Vec>>>, entity_types: Vec, + vertex_indices_to_ids: Vec, + vertex_ids_to_indices: HashMap, + cell_indices_to_ids: Vec, + cell_ids_to_indices: HashMap, } unsafe impl Sync for SerialSingleElementTopology {} impl SerialSingleElementTopology { - pub fn new(cells_input: &[usize], cell_type: ReferenceCellType) -> Self { + pub fn new( + cells_input: &[usize], + cell_type: ReferenceCellType, + point_indices_to_ids: &[usize], + grid_cell_indices_to_ids: &[usize], + ) -> Self { let size = reference_cell::entity_counts(cell_type)[0]; let ncells = cells_input.len() / size; + let mut vertex_indices_to_ids = vec![]; + let mut vertex_ids_to_indices = HashMap::new(); + let mut cell_indices_to_ids = vec![]; + let mut cell_ids_to_indices = HashMap::new(); + let mut index_map = vec![0; ncells]; let mut vertices = vec![]; let dim = reference_cell::dim(cell_type); @@ -59,11 +74,15 @@ impl SerialSingleElementTopology { for (cell_i, i) in index_map.iter_mut().enumerate() { let cell = &cells_input[start..start + size]; *i = cell_i; + cell_indices_to_ids.push(grid_cell_indices_to_ids[cell_i]); + cell_ids_to_indices.insert(grid_cell_indices_to_ids[cell_i], cell_i); let mut row = vec![]; for v in cell { if !vertices.contains(v) { entities_to_cells[0].push(vec![]); vertices.push(*v); + vertex_indices_to_ids.push(point_indices_to_ids[*v]); + vertex_ids_to_indices.insert(point_indices_to_ids[*v], *v); } row.push(vertices.iter().position(|&r| r == *v).unwrap()); } @@ -118,6 +137,10 @@ impl SerialSingleElementTopology { cells_to_entities, entities_to_cells, entity_types, + vertex_indices_to_ids, + vertex_ids_to_indices, + cell_indices_to_ids, + cell_ids_to_indices, } } } @@ -192,17 +215,17 @@ impl Topology for SerialSingleElementTopology { } } - fn vertex_index_to_id(&self, index: Self::IndexType) -> usize { - panic!(); + fn vertex_index_to_id(&self, index: usize) -> usize { + self.vertex_indices_to_ids[index] } - fn cell_index_to_id(&self, index: Self::IndexType) -> usize { - panic!(); + fn cell_index_to_id(&self, index: usize) -> usize { + self.cell_indices_to_ids[index] } - fn vertex_id_to_index(&self, id: usize) -> Self::IndexType { - panic!(); + fn vertex_id_to_index(&self, id: usize) -> usize { + self.vertex_ids_to_indices[&id] } - fn cell_id_to_index(&self, id: usize) -> Self::IndexType { - panic!(); + fn cell_id_to_index(&self, id: usize) -> usize { + self.cell_ids_to_indices[&id] } } @@ -211,7 +234,12 @@ mod test { use super::*; fn example_topology() -> SerialSingleElementTopology { - SerialSingleElementTopology::new(&[0, 1, 2, 2, 1, 3], ReferenceCellType::Triangle) + SerialSingleElementTopology::new( + &[0, 1, 2, 2, 1, 3], + ReferenceCellType::Triangle, + &[0, 1, 2, 3], + &[0, 1], + ) } #[test] diff --git a/src/grid/traits_impl.rs b/src/grid/traits_impl.rs index a640667..ca7dc5a 100644 --- a/src/grid/traits_impl.rs +++ b/src/grid/traits_impl.rs @@ -84,9 +84,7 @@ where type Geometry<'a> = CellGeometry<'a, T, GridImpl> where Self: 'a; fn id(&self) -> usize { - panic!(); - - //self.grid.geometry().point_index_to_id(self.index) + self.grid.geometry().point_index_to_id(self.index) } fn index(&self) -> usize { self.index