cmtool_data/
flowmap.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3use ndarray::Array2;
4
5use crate::{DataError, RawData, RawFlux, rawdata};
6
7pub struct FlowMapDescriptor {
8    pub flowmap: Array2<f64>,
9    pub neighbors: Array2<usize>, //TODO
10    pub volumes: Vec<f64>,
11}
12
13impl FlowMapDescriptor {
14    pub fn from_path(
15        data_flows_path: impl AsRef<std::path::Path>,
16        data_volumes_path: impl AsRef<std::path::Path>,
17    ) -> Result<Self, DataError> {
18        let df = rawdata::RawDataFlux::read_raw(data_flows_path).ok_or(DataError::BadData)?;
19        let dv = rawdata::RawDataScalar::read_raw(data_volumes_path).ok_or(DataError::BadData)?;
20
21        Self::from_raw_data(&df, &dv)
22    }
23
24    pub fn from_raw_data(
25        data_flows: &rawdata::RawDataFlux,
26        data_volumes: &rawdata::RawDataScalar,
27    ) -> Result<Self, DataError> {
28        let n_zone = data_flows.header.n_zone as usize;
29        let mut flowmap = Array2::<f64>::zeros((n_zone, n_zone));
30
31        let mut neighbors: Vec<Vec<usize>> = vec![Vec::with_capacity(10); n_zone];
32
33        if data_flows.header.n_zone != data_volumes.header.n_zone {
34            return Err(DataError::BadData);
35        }
36
37        let mut add_at = |i: usize, j: usize, val: f64| -> Result<(), DataError> {
38            //Error should never be triggered, by construction i<n and j<n
39            let g = flowmap.get_mut((i, j)).ok_or(DataError::BadData)?;
40            *g += val;
41            Ok(())
42        };
43
44        for &RawFlux {
45            id_source,
46            id_target,
47            flux_source_target,
48            flux_target_source,
49        } in data_flows.fluxes.iter()
50        {
51            let id_source = id_source as usize;
52            let id_target = id_target as usize;
53            assert!(flux_source_target >= 0.);
54            assert!(flux_target_source >= 0.);
55
56            add_at(id_source, id_target, flux_source_target)?;
57            add_at(id_target, id_source, flux_target_source)?;
58
59            neighbors[id_source].push(id_target);
60            neighbors[id_target].push(id_source);
61        }
62
63        let max_size = neighbors
64            .iter()
65            .map(|val| val.len())
66            .max()
67            .ok_or(DataError::BadData)?;
68
69        let mut neighbor_flat = Array2::<usize>::zeros((n_zone, max_size));
70
71        neighbor_flat.fill(n_zone + 1); //Any ghost neighbor will have value n+1
72
73        for (i_zone, neighbors_for_zone) in neighbors.iter().enumerate() {
74            for (i_n, id_neighbor) in neighbors_for_zone.iter().enumerate() {
75                // *(neighbor_flat.get_mut((i_zone, i_n)).unwrap()) = *id_neighbor;
76                *(neighbor_flat
77                    .get_mut((i_zone, i_n))
78                    .expect("Flat neighbor out of bound")) = *id_neighbor;
79            }
80        }
81
82        let volumes: Vec<f64> = data_volumes.values.iter().map(|v| v.value).collect();
83
84        Ok(FlowMapDescriptor {
85            flowmap,
86            neighbors: neighbor_flat,
87            volumes,
88        })
89    }
90}
91
92#[cfg(test)]
93mod test {
94    use super::*;
95    use crate::test_utils::chain_of_three;
96
97    fn check_shape(descriptor: &FlowMapDescriptor) {
98        assert!(!descriptor.volumes.is_empty());
99        assert!(descriptor.flowmap.is_square());
100        assert!(descriptor.volumes.len() == descriptor.flowmap.ncols());
101        assert!(descriptor.neighbors.nrows() == descriptor.flowmap.ncols());
102    }
103
104    #[test]
105    fn descriptor_from_synthetic_data() {
106        let (flow, volume) = chain_of_three();
107        let descriptor = FlowMapDescriptor::from_raw_data(&flow, &volume).unwrap();
108
109        check_shape(&descriptor);
110
111        // Flow map is indexed [from, to].
112        assert_eq!(descriptor.flowmap[[0, 1]], 1.0);
113        assert_eq!(descriptor.flowmap[[1, 0]], 0.5);
114        assert_eq!(descriptor.flowmap[[1, 2]], 2.0);
115        assert_eq!(descriptor.flowmap[[2, 1]], 0.25);
116        // No direct interface between the two ends of the chain.
117        assert_eq!(descriptor.flowmap[[0, 2]], 0.0);
118        assert_eq!(descriptor.flowmap[[2, 0]], 0.0);
119
120        assert_eq!(descriptor.volumes, vec![1.0, 2.0, 4.0]);
121
122        // Middle compartment has two neighbors, the ends have one plus a ghost.
123        let ghost = descriptor.volumes.len() + 1;
124        assert_eq!(descriptor.neighbors.ncols(), 2);
125        assert_eq!(descriptor.neighbors[[0, 0]], 1);
126        assert_eq!(descriptor.neighbors[[0, 1]], ghost);
127        assert_eq!(descriptor.neighbors[[1, 0]], 0);
128        assert_eq!(descriptor.neighbors[[1, 1]], 2);
129        assert_eq!(descriptor.neighbors[[2, 0]], 1);
130        assert_eq!(descriptor.neighbors[[2, 1]], ghost);
131    }
132
133    #[test]
134    fn descriptor_rejects_mismatched_zone_count() {
135        let (flow, _) = chain_of_three();
136        let short_volume = vec![1.0, 2.0].into();
137
138        assert!(FlowMapDescriptor::from_raw_data(&flow, &short_volume).is_err());
139    }
140
141    #[test]
142    fn read_descriptor() {
143        let _flow_cma = std::env::var("CUVE_SLDMSH_FLOW_PATH");
144
145        let _volume_cma = std::env::var("CUVE_SLDMSH_VOLUME_PATH");
146
147        if let (Ok(flow_cma), Ok(volume_cma)) = (_flow_cma, _volume_cma) {
148            let descriptor = FlowMapDescriptor::from_path(flow_cma, volume_cma).unwrap();
149
150            check_shape(&descriptor);
151        }
152    }
153}
cmtool_data/flowmap.rs

cmtool_data/
flowmap.rs
