Add implementation of grid topology and geometry

.github/workflows/style-checks.yml

@@ -23,7 +23,7 @@ jobs:
         uses: mpi4py/setup-mpi@v1
         with:
           mpi: mpich
-      -  uses: actions/checkout@v3
+      - uses: actions/checkout@v3
       - name: Install LAPACK & OpenBLAS
         run:
           sudo apt-get install libopenblas-dev liblapack-dev

Cargo.toml

@@ -6,9 +6,14 @@ edition = "2021"
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
+approx = "0.5"
+bempp-traits = { git = "https://github.com/bempp/bempp-rs.git" }
+bempp-tools = { git = "https://github.com/bempp/bempp-rs.git" }
+bempp-element = { git = "https://github.com/bempp/bempp-rs.git" }
+bempp-quadrature = { git = "https://github.com/bempp/bempp-rs.git" }
 num = "0.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"

src/grid.rs

@@ -1 +1,12 @@
-pub mod surface_triangle_grid;
+//! Grid implementation
+
+pub mod common;
+pub mod flat_triangle_grid;
+pub mod mixed_grid;
+pub mod single_element_grid;
+pub mod traits;
+pub mod traits_impl;
+
+pub use self::traits::Geometry;
+pub use self::traits::Grid;
+pub use self::traits::Topology;

src/grid/common.rs

@@ -0,0 +1,148 @@
+//! Functionality common to multiple grid implementations
+
+use crate::grid::traits::Geometry;
+use num::Float;
+use rlst_common::types::Scalar;
+use rlst_dense::{
+    array::Array,
+    traits::{Shape, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue},
+};
+
+/// Compute a physical point
+pub fn compute_point<T: Scalar, Table: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>>(
+    geometry: &impl Geometry<T = T>,
+    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;
+        }
+    }
+}
+
+/// Compute a Jacobian
+pub fn compute_jacobian<T: Scalar, Table: UnsafeRandomAccessByRef<4, Item = T> + Shape<4>>(
+    geometry: &impl Geometry<T = T>,
+    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]) };
+            }
+        }
+    }
+}
+
+/// Compute a normal from a Jacobian of a cell with topological dimension 2 and geometric dimension 3
+pub fn compute_normal_from_jacobian23<T: Scalar>(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::<T>());
+    for i in normal.iter_mut() {
+        *i /= size;
+    }
+}
+
+/// Compute the determinant of a matrix
+pub fn compute_det<T: Scalar<Real = T>>(jacobian: &[T], tdim: usize, gdim: usize) -> T {
+    assert_eq!(jacobian.len(), tdim * gdim);
+    match tdim {
+        1 => match gdim {
+            1 => T::abs(jacobian[0]),
+            2 => T::sqrt(jacobian.iter().map(|x| x.powi(2)).sum()),
+            3 => T::sqrt(jacobian.iter().map(|x| x.powi(2)).sum()),
+            _ => {
+                unimplemented!("compute_det() not implemented for topological dimension {tdim} and geometric dimension: {gdim}");
+            }
+        },
+        2 => match gdim {
+            2 => T::abs(jacobian[0] * jacobian[3] - jacobian[1] * jacobian[2]),
+            3 => T::sqrt(
+                [(1, 2), (2, 0), (0, 1)]
+                    .iter()
+                    .map(|(j, k)| {
+                        (jacobian[*j] * jacobian[3 + *k] - jacobian[*k] * jacobian[3 + *j]).powi(2)
+                    })
+                    .sum(),
+            ),
+            _ => {
+                unimplemented!("compute_det() not implemented for topological dimension {tdim} and geometric dimension: {gdim}");
+            }
+        },
+        3 => match gdim {
+            3 => T::abs(
+                [(0, 1, 2), (1, 2, 0), (2, 0, 1)]
+                    .iter()
+                    .map(|(i, j, k)| {
+                        jacobian[*i]
+                            * (jacobian[3 + *j] * jacobian[6 + *k]
+                                - jacobian[3 + *k] * jacobian[6 + *j])
+                    })
+                    .sum(),
+            ),
+            _ => {
+                unimplemented!("compute_det() not implemented for topological dimension {tdim} and geometric dimension: {gdim}");
+            }
+        },
+        _ => {
+            unimplemented!("compute_det() not implemented for topological dimension {tdim}");
+        }
+    }
+}
+
+/// Compute the diameter of a triangle
+pub fn compute_diameter_triangle<
+    T: Scalar<Real = T> + Float,
+    ArrayImpl: UnsafeRandomAccessByValue<1, Item = T> + Shape<1>,
+>(
+    v0: Array<T, ArrayImpl, 1>,
+    v1: Array<T, ArrayImpl, 1>,
+    v2: Array<T, ArrayImpl, 1>,
+) -> T {
+    let a = (v0.view() - v1.view()).norm_2();
+    let b = (v0 - v2.view()).norm_2();
+    let c = (v1 - v2).norm_2();
+    Scalar::sqrt((b + c - a) * (a + c - b) * (a + b - c) / (a + b + c))
+}
+
+/// Compute the diameter of a quadrilateral
+pub fn compute_diameter_quadrilateral<
+    T: Scalar<Real = T> + Float,
+    ArrayImpl: UnsafeRandomAccessByValue<1, Item = T> + Shape<1>,
+>(
+    v0: Array<T, ArrayImpl, 1>,
+    v1: Array<T, ArrayImpl, 1>,
+    v2: Array<T, ArrayImpl, 1>,
+    v3: Array<T, ArrayImpl, 1>,
+) -> T {
+    T::max((v0 - v3).norm_2(), (v1 - v2).norm_2())
+}

src/grid/flat_triangle_grid.rs

@@ -0,0 +1,9 @@
+//! A flat triangle grid
+//!
+//! In this grid, every cell is a flat triangle
+
+mod builder;
+mod grid;
+
+pub use self::builder::SerialFlatTriangleGridBuilder;
+pub use self::grid::SerialFlatTriangleGrid;

src/grid/flat_triangle_grid/builder.rs

@@ -0,0 +1,92 @@
+//! Grid builder
+
+use crate::grid::flat_triangle_grid::grid::SerialFlatTriangleGrid;
+use crate::traits::builder::Builder;
+use num::Float;
+use rlst_common::types::Scalar;
+use rlst_dense::{
+    array::{views::ArrayViewMut, Array},
+    base_array::BaseArray,
+    data_container::VectorContainer,
+    rlst_array_from_slice2, rlst_dynamic_array2,
+    traits::MatrixInverse,
+};
+use std::collections::HashMap;
+
+/// Grid builder for a flat triangle grid
+pub struct SerialFlatTriangleGridBuilder<T: Float + Scalar<Real = T>> {
+    points: Vec<T>,
+    cells: Vec<usize>,
+    point_indices_to_ids: Vec<usize>,
+    point_ids_to_indices: HashMap<usize, usize>,
+    cell_indices_to_ids: Vec<usize>,
+    cell_ids_to_indices: HashMap<usize, usize>,
+}
+
+impl<T: Float + Scalar<Real = T>> Builder<3> for SerialFlatTriangleGridBuilder<T>
+where
+    for<'a> Array<T, ArrayViewMut<'a, T, BaseArray<T, VectorContainer<T>, 2>, 2>, 2>: MatrixInverse,
+{
+    type GridType = SerialFlatTriangleGrid<T>;
+    type T = T;
+    type CellData = [usize; 3];
+    type GridMetadata = ();
+
+    fn new(_data: ()) -> Self {
+        Self {
+            points: vec![],
+            cells: vec![],
+            point_indices_to_ids: vec![],
+            point_ids_to_indices: HashMap::new(),
+            cell_indices_to_ids: vec![],
+            cell_ids_to_indices: HashMap::new(),
+        }
+    }
+
+    fn new_with_capacity(npoints: usize, ncells: usize, _data: ()) -> Self {
+        Self {
+            points: Vec::with_capacity(npoints * Self::GDIM),
+            cells: Vec::with_capacity(ncells * 3),
+            point_indices_to_ids: Vec::with_capacity(npoints),
+            point_ids_to_indices: HashMap::new(),
+            cell_indices_to_ids: Vec::with_capacity(ncells),
+            cell_ids_to_indices: HashMap::new(),
+        }
+    }
+
+    fn add_point(&mut self, id: usize, data: [T; 3]) {
+        self.point_ids_to_indices
+            .insert(id, self.point_indices_to_ids.len());
+        self.point_indices_to_ids.push(id);
+        self.points.extend_from_slice(&data);
+    }
+
+    fn add_cell(&mut self, id: usize, cell_data: [usize; 3]) {
+        self.cell_ids_to_indices
+            .insert(id, self.cell_indices_to_ids.len());
+        self.cell_indices_to_ids.push(id);
+        for id in &cell_data {
+            self.cells.push(self.point_ids_to_indices[id]);
+        }
+    }
+
+    fn create_grid(self) -> Self::GridType {
+        // TODO: remove this transposing
+        let npts = self.point_indices_to_ids.len();
+        let mut points = rlst_dynamic_array2!(T, [npts, 3]);
+        points.fill_from(rlst_array_from_slice2!(
+            T,
+            &self.points,
+            [npts, 3],
+            [1, npts]
+        ));
+        SerialFlatTriangleGrid::new(
+            points,
+            &self.cells,
+            self.point_indices_to_ids,
+            self.point_ids_to_indices,
+            self.cell_indices_to_ids,
+            self.cell_ids_to_indices,
+        )
+    }
+}