cmtool_core/model/
geometry.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3use std::{collections::BTreeSet, sync::Arc};
4
5use crate::{
6    coordinates::CartesianCoordinates,
7    ensight_gold::{self, types::ElementsType},
8    grid::{
9        CompartmentMesh, CylindricalAxis, MeshType, NeighborDirection, cylindrical_index, get_mesh,
10    },
11    model::{
12        CountVolumeElement,
13        data::{VerticesData, VolumeElementData},
14    },
15    utils::compute_centroid,
16};
17
18pub struct CMGeometry {
19    pub vertices: VerticesData,
20    pub volume_elements: VolumeElementData,
21    grid: Option<Box<dyn CompartmentMesh>>,
22    pub mesh_type: crate::grid::MeshType,
23}
24
25//Mutable
26impl CMGeometry {
27    fn fill_detail(&mut self, geometry: &Arc<ensight_gold::Geometry>) -> (Vec<usize>, Vec<usize>) {
28        let n_number_type = ensight_gold::types::VolumeElementTypes::NUMBER_OF_TYPES;
29        let n_part = geometry.number_of_part();
30        // let mut vertex_detail = Vec::<usize>::with_capacity(n_part);
31
32        let mut vertex_detail = vec![0; n_part];
33        let mut velem_detail = vec![0; n_part * n_number_type];
34        let mut n_vertex_total = 0;
35        let mut n_volume_elements_total = 0;
36
37        for (i_part, part) in geometry.parts.iter().enumerate() {
38            let base_index = i_part * n_number_type;
39            let part_n_vertex = part.n_vertex;
40            n_vertex_total += part_n_vertex;
41            vertex_detail[i_part] = part_n_vertex;
42            // vertex_detail.push(part.n_vertex);
43
44            for element in &part.elements {
45                if let ElementsType::VolumeElementType(vol_element) = element.etype {
46                    let index_element = vol_element.to_index();
47                    n_volume_elements_total += element.n_elements;
48                    //TODO:  prefetch slices in the outer loop as number of element type is known
49                    // Take mut slice [base_index.base_index+NUMBER_OF_TYPES]
50                    //Allow compiler to have linear indexing into slice: slice[index_element]+=elemenent.n_element
51                    velem_detail[base_index + index_element] += element.n_elements;
52                }
53            }
54        }
55
56        self.vertices.resize(n_vertex_total, &vertex_detail);
57
58        self.volume_elements
59            .resize(n_part, n_volume_elements_total, &velem_detail);
60
61        (vertex_detail, velem_detail)
62    }
63
64    fn init_cm_grid(&mut self, n_div: [usize; 3], mesh_type: MeshType) {
65        let mut axis: [crate::grid::AxisDescriptor; 3] = Default::default();
66
67        axis.iter_mut().zip(n_div).for_each(|(ax, div)| {
68            ax.n_range = div;
69        });
70
71        for vertex_global_id in 0..self.vertices.n_vertex() {
72            axis.iter_mut()
73                .zip(self.vertices.get_slice_xyz(vertex_global_id))
74                .for_each(|(axe, &vertex)| {
75                    axe.min_range = axe.min_range.min(vertex);
76                    axe.max_range = axe.max_range.max(vertex);
77                });
78            if mesh_type == MeshType::Cylindrical {
79                let radius = self.vertices.get_radius_from_global_id(vertex_global_id);
80                axis[cylindrical_index(CylindricalAxis::R)].max_range =
81                    axis[0].max_range.max(radius);
82            }
83        }
84        if mesh_type == MeshType::Cylindrical {
85            axis[cylindrical_index(CylindricalAxis::R)].min_range = 0.;
86        }
87
88        // axis.iter_mut()
89        //     .for_each(|ax| ax.step = (ax.max_range - ax.min_range) / (ax.n_range as f64));
90
91        self.grid = Some(get_mesh(mesh_type, axis));
92    }
93
94    fn detect_compartment(&mut self, n_div: [usize; 3], mesh_type: MeshType) {
95        self.init_cm_grid(n_div, mesh_type);
96
97        let grid = self.grid.as_ref().unwrap();
98
99        let vertices_id: Vec<_> = (0..self.vertices.n_vertex())
100            .map(|global_id| {
101                grid.cell_from_coordinates(self.vertices.get_slice_xyz(global_id))
102                    .unwrap()
103            })
104            .collect();
105
106        for vol_element_global_id in 0..self.volume_elements.n_element() {
107            let n_vertex = self
108                .volume_elements
109                .get_vertex_per_element(vol_element_global_id);
110
111            self.volume_elements.xyz[vol_element_global_id] =
112                self.get_element_centroid(vol_element_global_id, n_vertex);
113
114            let unique_cids: BTreeSet<_> = (0..n_vertex)
115                .map(|k_vertex| {
116                    let vertex_global_id = self
117                        .volume_elements
118                        .get_vertex_from_vol_global_id(vol_element_global_id, k_vertex);
119                    // self.vertices.ve_id[vertex_global_id]
120                    vertices_id[vertex_global_id]
121                })
122                .collect();
123
124            self.volume_elements
125                .set_number_cid(vol_element_global_id, unique_cids.len());
126
127            for (index, c_id) in unique_cids.into_iter().enumerate() {
128                self.volume_elements
129                    .set_list_compartment_id(vol_element_global_id, index, c_id);
130            }
131        }
132    }
133
134    pub fn get_element_centroid(
135        &self,
136        vol_element_global_id: usize,
137        n_vertex: usize,
138    ) -> CartesianCoordinates {
139        //This iter is lazy as we map without operation
140        let iter = (0..n_vertex).map(|k| {
141            let vertex_id = self
142                .volume_elements
143                .get_vertex_from_vol_global_id(vol_element_global_id, k);
144            self.vertices.get_slice_xyz(vertex_id)
145        });
146
147        compute_centroid(iter)
148    }
149}
150
151impl CMGeometry {
152    pub fn n_zone(&self) -> usize {
153        self.grid.as_ref().unwrap().number_cell()
154    }
155
156    pub(super) fn interface_iter(&self) -> impl Iterator<Item = (usize, usize, usize, usize)> + '_ {
157        (0..self.volume_elements.n_element()).flat_map(move |vol_element_global_id| {
158            let interface_cid_0 = self
159                .volume_elements
160                .get_list_compartment_id(vol_element_global_id, 0);
161
162            let n_cid = self.volume_elements.get_number_cid(vol_element_global_id);
163
164            (0..n_cid).map(move |k_vertex| {
165                let interface_cid_k = self
166                    .volume_elements
167                    .get_list_compartment_id(vol_element_global_id, k_vertex);
168
169                (
170                    vol_element_global_id,
171                    interface_cid_0,
172                    interface_cid_k,
173                    k_vertex,
174                )
175            })
176        })
177    }
178
179    pub fn fill_vertices(
180        &self,
181        volume_element_global_id: usize,
182        n_vertex: usize,
183        local_vertices: &mut Vec<CartesianCoordinates>,
184    ) {
185        local_vertices.clear();
186        local_vertices.reserve(n_vertex);
187
188        for k_vertex in 0..n_vertex {
189            let vertex_global_id = self
190                .volume_elements
191                .get_vertex_from_vol_global_id(volume_element_global_id, k_vertex);
192
193            local_vertices.push(CartesianCoordinates(
194                self.vertices.get_slice_xyz(vertex_global_id).to_owned(),
195            ));
196        }
197    }
198
199    pub fn get_count_volume_element_first_pass(&self) -> CountVolumeElement {
200        let mut count = CountVolumeElement::new(self.n_zone());
201
202        let grid = self.grid.as_ref().unwrap();
203
204        for (vol_element_global_id, _, interface_cid_k, k_vertex) in self.interface_iter() {
205            count.incr_compartment(interface_cid_k);
206            // if k_vertex >= 1 {
207            for i in 0..k_vertex {
208                let cid_i = self
209                    .volume_elements
210                    .get_list_compartment_id(vol_element_global_id, i);
211                let neighbors = grid.are_cell_neighbor(cid_i, interface_cid_k);
212                if neighbors != NeighborDirection::NotNeighbors {
213                    let (id1, id2) = neighbors.ordered_pair(cid_i, interface_cid_k);
214                    count.incr_interface(id1, id2);
215                }
216            }
217            // }
218        }
219
220        count
221    }
222
223    pub fn init(
224        n_div: [usize; 3],
225        geometry: Arc<ensight_gold::Geometry>,
226        mesh_type: crate::grid::MeshType,
227    ) -> Self {
228        let mut cm_geometry = Self {
229            vertices: Default::default(),
230            volume_elements: Default::default(),
231            grid: None,
232            mesh_type,
233        };
234
235        let (vertex_detail, velem_detail) = cm_geometry.fill_detail(&geometry);
236
237        let mut global_vertex_counter = 0;
238        let mut global_volume_element_counter = 0;
239
240        for part_it in geometry.parts.iter().enumerate() {
241            global_vertex_counter +=
242                cm_geometry
243                    .vertices
244                    .fill_from_part(global_vertex_counter, &vertex_detail, part_it);
245
246            global_volume_element_counter += cm_geometry.volume_elements.fill_from_part(
247                global_volume_element_counter,
248                part_it,
249                &velem_detail,
250                &cm_geometry.vertices,
251            );
252        }
253
254        cm_geometry.detect_compartment(n_div, mesh_type);
255
256        cm_geometry
257    }
258
259    pub fn get_grid(&self) -> Option<&dyn CompartmentMesh> {
260        self.grid.as_deref()
261    }
262}
cmtool_core/model/geometry.rs

cmtool_core/model/
geometry.rs
