update weekly workflow

.github/workflows/run-weekly-tests.yml

@@ -3,6 +3,8 @@ name: 🧪
 on:
   schedule:
     - cron: "0 7 * * 1"
+  push:
+    branches: ["**"]
 
 jobs:
   run-tests-rust:
@@ -12,6 +14,9 @@ jobs:
       matrix:
         rust-version: ["stable", "beta", "nightly"]
         mpi: [ 'mpich', 'openmpi']
+        f_mpi: ["", "mpi,"]
+        f_strict: ["", "strict,"]
+        f_serde: ["", "serde,"]
     steps:
       - name: Set up Rust
         uses: actions-rust-lang/setup-rust-toolchain@v1
@@ -30,30 +35,23 @@ jobs:
           sudo apt-get install -y libopenblas-dev liblapack-dev
 
       - name: Build rust library
-        run: cargo build --features "strict"
+        run: cargo build --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
       - name: Build rust library in release mode
-        run: cargo build --release --features "strict"
-      - name: Build rust library in release mode with mpi
-        run: cargo build --release --features "strict,mpi"
-
+        run: cargo build --release --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
 
       - name: Run unit tests
-        run: cargo test --features "strict"
+        run: cargo test --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
       - name: Run unit tests in release mode
-        run: cargo test --release --features "strict"
-      - name: Run unit tests with mpi enabled
-        run: cargo test --features "mpi,strict"
-      - name: Run unit tests in release mode with mpi enabled
-        run: cargo test --release --features "mpi,strict"
+        run: cargo test --release --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
       - name: Run tests
-        run: cargo test --examples --release --features "mpi,strict"
-      # - name: Test benchmarks build
-      #   run: cargo bench --no-run --features "mpi,strict"
+        run: cargo test --examples --release --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
+      - name: Test benchmarks build
+        run: cargo bench --no-run --features "${{f_mpi}}${{f_strict}}${{f_serde}}"
       - name: Run examples
         run: |
           python3 find_examples.py
           chmod +x examples.sh
           ./examples.sh
 
       - name: Build docs
-        run: cargo doc --features "mpi,strict" --no-deps
+        run: cargo doc --features "${{f_mpi}}${{f_strict}}${{f_serde}}" --no-deps

.gitignore

@@ -13,7 +13,11 @@ Cargo.lock
 
 examples.sh
 
-_*.msh
 *.btr
 
 .vscode/
+
+# Files generated by tests
+_*.msh
+_*.ron
+

Cargo.toml

@@ -1,6 +1,6 @@
 [features]
 mpi = ["dep:mpi", "dep:coupe", "ndelement/mpi"]
-serde = ["dep:serde", "ndelement/serde"]
+serde = ["dep:serde", "ndelement/serde", "dep:ron"]
 strict = []
 
 [package]
@@ -28,6 +28,7 @@ ndelement = "0.1.0"
 num = "0.4"
 rlst = "0.2.0"
 serde = { version = "1", features = ["derive"], optional = true }
+ron = { version = "0.8", optional = true }
 
 [dev-dependencies]
 approx = "0.5"

examples/test_parallel_io.rs

@@ -0,0 +1,99 @@
+//? mpirun -n {{NPROCESSES}} --features "mpi,serde"
+
+#[cfg(feature = "mpi,serde")]
+use mpi::{collective::CommunicatorCollectives, environment::Universe, traits::Communicator};
+#[cfg(feature = "mpi,serde")]
+use ndelement::{ciarlet::CiarletElement, types::ReferenceCellType};
+#[cfg(feature = "mpi,serde")]
+use ndgrid::{
+    grid::parallel::ParallelGrid,
+    traits::{Builder, Grid, ParallelBuilder, RONExportParallel, RONImportParallel},
+    SingleElementGrid, SingleElementGridBuilder,
+};
+
+#[cfg(feature = "mpi,serde")]
+fn create_single_element_grid_data(b: &mut SingleElementGridBuilder<f64>, n: usize) {
+    for y in 0..n {
+        for x in 0..n {
+            b.add_point(
+                y * n + x,
+                &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64, 0.0],
+            );
+        }
+    }
+
+    for i in 0..n - 1 {
+        for j in 0..n - 1 {
+            b.add_cell(
+                i * (n - 1) + j,
+                &[j * n + i, j * n + i + 1, j * n + i + n, j * n + i + n + 1],
+            );
+        }
+    }
+}
+
+#[cfg(feature = "mpi,serde")]
+fn example_single_element_grid<C: Communicator>(
+    comm: &C,
+    n: usize,
+) -> ParallelGrid<'_, C, SingleElementGrid<f64, CiarletElement<f64>>> {
+    let rank = comm.rank();
+
+    let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
+
+    if rank == 0 {
+        create_single_element_grid_data(&mut b, n);
+        b.create_parallel_grid(comm)
+    } else {
+        b.receive_parallel_grid(comm, 0)
+    }
+}
+
+#[cfg(feature = "mpi,serde")]
+fn test_parallel_export<C: Communicator>(comm: &C) {
+    let size = comm.size();
+
+    let n = 10;
+    let grid = example_single_element_grid(comm, n);
+    let filename = format!("_examples_parallel_io_{}ranks.ron", size);
+    grid.export_as_ron(&filename);
+}
+
+#[cfg(feature = "mpi,serde")]
+fn test_parallel_import<C: Communicator>(comm: &C) {
+    let size = comm.size();
+
+    let filename = format!("_examples_parallel_io_{}ranks.ron", size);
+    let grid = ParallelGrid::<'_, C, SingleElementGrid<f64, CiarletElement<f64>>>::import_from_ron(
+        comm, &filename,
+    );
+
+    let n = 10;
+    let grid2 = example_single_element_grid(comm, n);
+
+    assert_eq!(
+        grid.entity_count(ReferenceCellType::Point),
+        grid2.entity_count(ReferenceCellType::Point)
+    );
+}
+
+#[cfg(feature = "mpi,serde")]
+fn main() {
+    let universe: Universe = mpi::initialize().unwrap();
+    let world = universe.world();
+    let rank = world.rank();
+
+    if rank == 0 {
+        println!("Testing parallel grid export");
+    }
+    test_parallel_export(&world);
+
+    world.barrier();
+
+    if rank == 0 {
+        println!("Testing parallel grid import");
+    }
+    test_parallel_import(&world);
+}
+#[cfg(not(feature = "mpi,serde"))]
+fn main() {}

src/geometry.rs

@@ -1,7 +1,7 @@
 //! Grid geometry
 mod geometry_map;
 mod point;
-mod single_element;
+pub(crate) mod single_element;
 pub use geometry_map::GeometryMap;
 pub use point::{Point, PointIter};
 pub use single_element::{SingleElementEntityGeometry, SingleElementGeometry};

src/geometry/geometry_map.rs

@@ -4,9 +4,8 @@ use crate::{
     types::{Array2D, ArrayND, RealScalar},
 };
 use ndelement::{reference_cell, traits::FiniteElement, types::ReferenceCellType};
-use rlst::{
-    rlst_dynamic_array2, rlst_dynamic_array4, DefaultIterator, RandomAccessByRef, RlstScalar, Shape,
-};
+use rlst::UnsafeRandomAccessByRef;
+use rlst::{rlst_dynamic_array4, RandomAccessByRef, RlstScalar, Shape};
 
 /// Single element geometry
 #[derive(Debug)]
@@ -18,26 +17,37 @@ pub struct GeometryMap<'a, T: RealScalar> {
     table: ArrayND<4, T>,
 }
 
-fn det<T: RlstScalar>(m: &Array2D<T>) -> T {
-    match m.shape()[0] {
-        0 => T::from(1.0).unwrap(),
-        1 => m[[0, 0]],
-        2 => m[[0, 0]] * m[[1, 1]] - m[[0, 1]] * m[[1, 0]],
-        3 => {
-            m[[0, 0]] * (m[[1, 1]] * m[[2, 2]] - m[[1, 2]] * m[[2, 1]])
-                - m[[0, 1]] * (m[[1, 0]] * m[[2, 2]] - m[[1, 2]] * m[[2, 0]])
-                + m[[0, 2]] * (m[[1, 0]] * m[[2, 1]] - m[[1, 1]] * m[[2, 0]])
+fn norm<T: RlstScalar>(vector: &[T]) -> T {
+    vector.iter().map(|&i| i * i).sum::<T>().sqrt()
+}
+
+fn cross<T: RlstScalar>(mat: &[T], result: &mut [T]) {
+    match mat.len() {
+        0 => {}
+        2 => {
+            debug_assert!(result.len() == 2);
+            unsafe {
+                *result.get_unchecked_mut(0) = *mat.get_unchecked(1);
+                *result.get_unchecked_mut(1) = -*mat.get_unchecked(0);
+            }
+        }
+        6 => {
+            debug_assert!(result.len() == 3);
+            unsafe {
+                *result.get_unchecked_mut(0) = *mat.get_unchecked(1) * *mat.get_unchecked(5)
+                    - *mat.get_unchecked(2) * *mat.get_unchecked(4);
+                *result.get_unchecked_mut(1) = *mat.get_unchecked(2) * *mat.get_unchecked(3)
+                    - *mat.get_unchecked(0) * *mat.get_unchecked(5);
+                *result.get_unchecked_mut(2) = *mat.get_unchecked(0) * *mat.get_unchecked(4)
+                    - *mat.get_unchecked(1) * *mat.get_unchecked(3);
+            }
         }
         _ => {
             unimplemented!();
         }
     }
 }
 
-fn norm<T: RlstScalar>(vector: &[T]) -> T {
-    vector.iter().map(|&i| i * i).sum::<T>().sqrt()
-}
-
 impl<'a, T: RealScalar> GeometryMap<'a, T> {
     /// Create new
     pub fn new<A2D: RandomAccessByRef<2, Item = T> + Shape<2>>(
@@ -47,7 +57,7 @@ impl<'a, T: RealScalar> GeometryMap<'a, T> {
         entities: &'a Array2D<usize>,
     ) -> Self {
         let tdim = reference_cell::dim(element.cell_type());
-        assert_eq!(points.shape()[0], tdim);
+        debug_assert!(points.shape()[0] == tdim);
         let gdim = geometry_points.shape()[0];
         let npoints = points.shape()[1];
 
@@ -78,35 +88,38 @@ impl<'a, T: RealScalar> GeometryMapTrait for GeometryMap<'a, T> {
     }
     fn points(&self, entity_index: usize, points: &mut [T]) {
         let npts = self.table.shape()[1];
-        assert_eq!(points.len(), self.gdim * npts);
+        debug_assert!(points.len() == self.gdim * npts);
 
-        points.fill(T::from(0.0).unwrap());
-        // TODO: can I use rlst better here?
+        points.fill(T::zero());
         for i in 0..self.entities.shape()[0] {
-            let v = self.entities[[i, entity_index]];
+            let v = unsafe { *self.entities.get_unchecked([i, entity_index]) };
             for point_index in 0..npts {
-                let t = *self.table.get([0, point_index, i, 0]).unwrap();
+                let t = unsafe { *self.table.get_unchecked([0, point_index, i, 0]) };
                 for gd in 0..self.gdim {
-                    points[gd + self.gdim * point_index] +=
-                        *self.geometry_points.get([gd, v]).unwrap() * t;
+                    unsafe {
+                        *points.get_unchecked_mut(gd + self.gdim * point_index) +=
+                            *self.geometry_points.get_unchecked([gd, v]) * t
+                    };
                 }
             }
         }
     }
     fn jacobians(&self, entity_index: usize, jacobians: &mut [T]) {
         let npts = self.table.shape()[1];
-        assert_eq!(jacobians.len(), self.gdim * self.tdim * npts);
+        debug_assert!(jacobians.len() == self.gdim * self.tdim * npts);
 
-        jacobians.fill(T::from(0.0).unwrap());
-        // TODO: can I use rlst better here?
+        jacobians.fill(T::zero());
         for i in 0..self.entities.shape()[0] {
-            let v = self.entities[[i, entity_index]];
+            let v = unsafe { *self.entities.get_unchecked([i, entity_index]) };
             for point_index in 0..npts {
-                for gd in 0..self.gdim {
-                    for td in 0..self.tdim {
-                        jacobians[gd + self.gdim * td + self.gdim * self.tdim * point_index] +=
-                            *self.geometry_points.get([gd, v]).unwrap()
-                                * *self.table.get([1 + td, point_index, i, 0]).unwrap();
+                for td in 0..self.tdim {
+                    let t = unsafe { *self.table.get_unchecked([1 + td, point_index, i, 0]) };
+                    for gd in 0..self.gdim {
+                        unsafe {
+                            *jacobians.get_unchecked_mut(
+                                gd + self.gdim * td + self.gdim * self.tdim * point_index,
+                            ) += *self.geometry_points.get_unchecked([gd, v]) * t
+                        };
                     }
                 }
             }
@@ -124,49 +137,27 @@ impl<'a, T: RealScalar> GeometryMapTrait for GeometryMap<'a, T> {
             panic!("Can only compute normal for entities where tdim + 1 == gdim");
         }
         let npts = self.table.shape()[1];
-        assert_eq!(jacobians.len(), self.gdim * self.tdim * npts);
-        assert_eq!(jdets.len(), npts);
-        assert_eq!(normals.len(), self.gdim * npts);
+        debug_assert!(jacobians.len() == self.gdim * self.tdim * npts);
+        debug_assert!(jdets.len() == npts);
+        debug_assert!(normals.len() == self.gdim * npts);
 
-        // TODO: can we remove this memory assignment?
-        let mut temp = rlst_dynamic_array2!(T, [self.tdim, self.tdim]);
+        self.jacobians(entity_index, jacobians);
 
         for point_index in 0..npts {
-            for gd in 0..self.gdim {
-                for td in 0..self.tdim {
-                    jacobians[gd + self.gdim * td + self.gdim * self.tdim * point_index] = self
-                        .entities
-                        .view()
-                        .slice(1, entity_index)
-                        .iter()
-                        .enumerate()
-                        .map(|(i, v)| {
-                            *self.geometry_points.get([gd, v]).unwrap()
-                                * *self.table.get([1 + td, point_index, i, 0]).unwrap()
-                        })
-                        .sum::<T>();
-                }
-            }
-            for gd in 1..self.gdim {
-                for td in 0..self.tdim {
-                    temp[[gd - 1, td]] =
-                        jacobians[gd + self.gdim * td + self.gdim * self.tdim * point_index];
-                }
-            }
-            for gd in 0..self.gdim {
-                normals[gd + self.gdim * point_index] =
-                    if gd % 2 == 0 { det(&temp) } else { -det(&temp) };
-                if gd < self.tdim {
-                    for td in 0..self.tdim {
-                        temp[[gd, td]] =
-                            jacobians[gd + self.gdim * td + self.gdim * self.tdim * point_index];
-                    }
-                }
-            }
-            jdets[point_index] =
-                norm(&normals[self.gdim * point_index..self.gdim * (point_index + 1)]);
-            for gd in 0..self.gdim {
-                normals[gd + self.gdim * point_index] /= jdets[point_index];
+            let j = unsafe {
+                &jacobians.get_unchecked(
+                    self.gdim * self.tdim * point_index..self.gdim * self.tdim * (point_index + 1),
+                )
+            };
+            let n = unsafe {
+                &mut normals
+                    .get_unchecked_mut(self.gdim * point_index..self.gdim * (point_index + 1))
+            };
+            let jd = unsafe { jdets.get_unchecked_mut(point_index) };
+            cross(j, n);
+            *jd = norm(n);
+            for n_i in n.iter_mut() {
+                *n_i /= *jd;
             }
         }
     }

src/geometry/single_element.rs

@@ -3,6 +3,8 @@ mod entity_geometry;
 mod geometry;
 
 pub use entity_geometry::SingleElementEntityGeometry;
+#[cfg(feature = "serde")]
+pub(crate) use geometry::SerializableGeometry;
 pub use geometry::SingleElementGeometry;
 
 #[cfg(test)]