cmtool_core/
lib.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3pub mod coordinates;
4pub mod ensight_gold;
5mod errors;
6pub mod grid;
7pub mod model;
8pub mod utils;
9pub use errors::CoreError;
10
11#[cfg(feature = "use_vtk")]
12use grid::vtk::VtkCm;
13#[cfg(feature = "use_vtk")]
14use grid::vtk::add_celldata_to_vtk;
15
16use cmtool_data::{RawData, RawDataFlux, RawDataScalar};
17use model::{CMGeometry, CMModel, Scalar, Vector};
18use std::{path::Path, sync::Arc};
19
20fn resolve_path(
21    root: &impl AsRef<std::path::Path>,
22    relative_path: &str,
23) -> impl AsRef<std::path::Path> {
24    std::path::PathBuf::from(root.as_ref()).join(relative_path)
25}
26
27pub enum ExportType {
28    EnsightGold,
29}
30
31// trait GeometryInfo {}
32
33// trait CfdCase {
34//     fn get_root(&self) -> String;
35//     fn get_geometry_relative_path(&self) -> String;
36// }
37
38pub struct CMHandle {
39    model: Arc<model::CMModel>,
40    ///Knobs of the balancing pass, the caller may replace them before generating
41    balance: model::BalanceSettings,
42    _root_result: String, //TODO EITHER USE IT OR REMOVE
43    eg_geometry: Arc<ensight_gold::Geometry>,
44    cm_geometry: Arc<CMGeometry>,
45}
46
47impl CMHandle {
48    ///Replaces the balancing knobs used when a flow map is generated
49    pub fn set_balance_settings(&mut self, settings: model::BalanceSettings) {
50        self.balance = settings;
51    }
52
53    pub fn balance_settings(&self) -> &model::BalanceSettings {
54        &self.balance
55    }
56
57    pub fn grid(&self) -> &dyn crate::grid::CompartmentMesh {
58        self.model.grid()
59    }
60
61    pub fn init(
62        n_div: [usize; 3],
63        root: &str,
64        geometry_filename: &str,
65        _meshtype: grid::MeshType,
66    ) -> Result<Self, CoreError> {
67        let fullpath = format!("{}/{}", root, geometry_filename);
68
69        let eg_geometry = Arc::new(ensight_gold::Geometry::new(Path::new(&fullpath))?);
70
71        println!("{}", eg_geometry);
72
73        let cm_geometry = Arc::new(CMGeometry::init(
74            n_div,
75            eg_geometry.clone(),
76            grid::MeshType::Cylindrical,
77        ));
78
79        Ok(Self {
80            model: Arc::new(CMModel::init(cm_geometry.clone())),
81            balance: Default::default(),
82            _root_result: String::from("./test"),
83            eg_geometry,
84            cm_geometry,
85        })
86    }
87
88    pub fn dump_volume(&self) {
89        self.model.compartments_volumes();
90
91        todo!()
92    }
93
94    pub fn dump_all(
95        &self,
96        root_export: impl AsRef<std::path::Path>,
97        root_input: impl AsRef<std::path::Path>,
98        vars: &[ensight_gold::case::VariableInfo],
99    ) -> Result<(), CoreError> {
100        std::fs::create_dir_all(&root_export)?;
101
102        //Scalars and vectors are dumped independently, only the scalars are kept for the vtk export
103        let mut rs = Vec::with_capacity(vars.len());
104
105        for v in vars.iter() {
106            match v.get_type() {
107                ensight_gold::case::VariableType::Scalar => {
108                    let sc = self.dump_scalar(
109                        resolve_path(&root_export, &v.name),
110                        resolve_path(&root_input, &v.filepath),
111                    )?;
112                    rs.push((sc, v.name.clone()));
113                }
114                ensight_gold::case::VariableType::Vector => {
115                    self.dump_vector(
116                        resolve_path(&root_export, &v.name),
117                        resolve_path(&root_input, &v.filepath),
118                    )?;
119                }
120            }
121        }
122
123        #[cfg(feature = "use_vtk")]
124        self.export_vtk(
125            format!("{}/cma_case.vtu", root_export.as_ref().display()),
126            rs,
127        )?;
128
129        Ok(())
130    }
131
132    pub fn get_scalar(&self, path: impl AsRef<std::path::Path>) -> Result<Scalar, CoreError> {
133        Self::s_get_scalar(path, self.eg_geometry.clone(), self.cm_geometry.clone())
134    }
135
136    fn s_get_scalar(
137        path: impl AsRef<std::path::Path>,
138        eg_geometry: Arc<ensight_gold::Geometry>,
139        cm_geometry: Arc<CMGeometry>,
140    ) -> Result<Scalar, CoreError> {
141        let s = ensight_gold::scalar::ScalarField::init(eg_geometry.clone(), path)?;
142        Ok(Scalar::new(s, &cm_geometry, &eg_geometry))
143    }
144
145    fn ts_dump_scalar(
146        res_name: impl AsRef<std::path::Path>,
147        path: impl AsRef<std::path::Path>,
148        eg_geometry: Arc<ensight_gold::Geometry>,
149        cm_geometry: Arc<CMGeometry>,
150        model: Arc<model::CMModel>,
151    ) -> Result<RawDataScalar, CoreError> {
152        let scalar = Self::s_get_scalar(path, eg_geometry, cm_geometry)?;
153
154        let scalar_data = model.export_volume_integral_per_zone(scalar)?;
155
156        scalar_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
157
158        Ok(scalar_data)
159    }
160
161    pub fn dump_scalar(
162        &self,
163        res_name: impl AsRef<std::path::Path>,
164        path: impl AsRef<std::path::Path>,
165    ) -> Result<RawDataScalar, CoreError> {
166        Self::ts_dump_scalar(
167            res_name,
168            path,
169            self.eg_geometry.clone(),
170            self.cm_geometry.clone(),
171            self.model.clone(),
172        )
173    }
174
175    ///Integrates a scalar over each compartment, weighting every element by `phase_fraction`
176    ///first, so a field carried by one phase is integrated over the volume that phase occupies.
177    ///Scaling the compartment integral afterwards is not the same number unless the field and
178    ///the fraction are uncorrelated inside the compartment.
179    pub fn dump_scalar_fraction(
180        &self,
181        res_name: impl AsRef<std::path::Path>,
182        path: impl AsRef<std::path::Path>,
183        phase_fraction: Scalar,
184    ) -> Result<RawDataScalar, CoreError> {
185        let scalar = self.get_scalar(path)?.element_wise(&phase_fraction)?;
186        let scalar_data = self.model.export_volume_integral_per_zone(scalar)?;
187
188        scalar_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
189
190        Ok(scalar_data)
191    }
192
193    pub fn dump_real_volume(&self, res_name: impl AsRef<std::path::Path>) -> Result<(), CoreError> {
194        let volumes_data: RawDataScalar = self.model.get_real_volume().into();
195
196        volumes_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
197
198        Ok(())
199    }
200
201    pub fn dump_vector_phase_fraction(
202        &self,
203        res_name: impl AsRef<std::path::Path>,
204        path: impl AsRef<std::path::Path>,
205        phase_fraction: Scalar,
206    ) -> Result<(), CoreError> {
207        let v = ensight_gold::vectors::VectorField::init(self.eg_geometry.clone(), path)?;
208        let vector =
209            Vector::new(v, &self.cm_geometry, &self.eg_geometry).scale_by(phase_fraction)?;
210
211        let flow_data = self
212            .model
213            .compute_flux_between_compartments(vector, &self.balance)?;
214
215        flow_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
216        Ok(())
217    }
218
219    pub fn dump_vector_raw(
220        &self,
221        res_name: impl AsRef<std::path::Path>,
222        vector: Vector,
223    ) -> Result<(), CoreError> {
224        let flow_data = self
225            .model
226            .compute_flux_between_compartments(vector, &self.balance)?;
227
228        flow_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
229        Ok(())
230    }
231
232    pub fn dump_vector(
233        &self,
234        res_name: impl AsRef<std::path::Path>,
235        path: impl AsRef<std::path::Path>,
236    ) -> Result<(), CoreError> {
237        let v = ensight_gold::vectors::VectorField::init(self.eg_geometry.clone(), path)?;
238        let vector = Vector::new(v, &self.cm_geometry, &self.eg_geometry);
239        let flow_data = self
240            .model
241            .compute_flux_between_compartments(vector, &self.balance)?;
242
243        flow_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
244        Ok(())
245    }
246
247    pub fn vector_from_scalar(
248        &self,
249        path_i: impl AsRef<std::path::Path>,
250        path_j: impl AsRef<std::path::Path>,
251        path_k: impl AsRef<std::path::Path>,
252    ) -> Result<Vector, CoreError> {
253        let s = ensight_gold::scalar::ScalarField::init(self.eg_geometry.clone(), path_i)?;
254        let sj = ensight_gold::scalar::ScalarField::init(self.eg_geometry.clone(), path_j)?;
255        let sk = ensight_gold::scalar::ScalarField::init(self.eg_geometry.clone(), path_k)?;
256        Vector::from_scalar([s, sj, sk], &self.cm_geometry, &self.eg_geometry)
257    }
258
259    ///Volume of each compartment, as covered by the mesh
260    pub fn real_volume(&self) -> Vec<f64> {
261        self.model.get_real_volume()
262    }
263
264    ///Flow map out of the three components of a velocity, optionally scaled by the volume
265    ///fraction of its phase: a phase only carries its own share of the flow
266    pub fn dump_vector_from_scalar(
267        &self,
268        res_name: impl AsRef<std::path::Path>,
269        path_i: impl AsRef<std::path::Path>,
270        path_j: impl AsRef<std::path::Path>,
271        path_k: impl AsRef<std::path::Path>,
272        phase_fraction: Option<Scalar>,
273    ) -> Result<RawDataFlux, CoreError> {
274        let vector = self.vector_from_scalar(path_i, path_j, path_k)?;
275        let vector = match phase_fraction {
276            Some(fraction) => vector.scale_by(fraction)?,
277            None => vector,
278        };
279        let flow_data = self
280            .model
281            .compute_flux_between_compartments(vector, &self.balance)?;
282
283        flow_data.write_raw(&format!("{}.raw", res_name.as_ref().to_str().unwrap()))?;
284        Ok(flow_data)
285    }
286
287    pub fn export_geometry_compartments(&self) {
288        todo!()
289    }
290
291    // #[cfg(feature = "use_vtk")]
292    // pub fn write_vtk(&self, path: impl AsRef<std::path::Path>) -> Result<(), CoreError> {
293    //     let mesh = self.cm_geometry.get_grid().unwrap();
294    //     let p = path.as_ref().to_str().unwrap();
295    //     let mut vtk = mesh.get_vtk(p)?;
296
297    //     let volumes_data = self.model.get_real_volume();
298
299    //     let volumes_data_array = vtkio::model::DataArray::scalars("real_volume", 1);
300
301    //     let volumes_data_array = volumes_data_array.with_vec(volumes_data);
302
303    //     add_celldata_to_vtk(
304    //         &mut vtk,
305    //         vtkio::model::Attribute::DataArray(volumes_data_array),
306    //     );
307
308    //     let mut vtk_bytes = Vec::<u8>::new();
309    //     vtk.write_xml(&mut vtk_bytes).unwrap();
310    //     std::fs::write(path, vtk_bytes).unwrap();
311
312    //     Ok(())
313    // }
314
315    #[cfg(feature = "use_vtk")]
316    pub fn export_vtk(
317        &self,
318        path: impl AsRef<std::path::Path>,
319        sc: Vec<(RawDataScalar, String)>,
320    ) -> Result<(), CoreError> {
321        let mesh = self.cm_geometry.get_grid().unwrap();
322        let p = path.as_ref().to_str().unwrap();
323        let mut vtk = mesh.get_vtk(p)?;
324
325        let mut add_cell = |name: &str, data: Vec<f64>| {
326            let data_array = vtkio::model::DataArray::scalars(name, 1);
327            let data_array = data_array.with_vec(data);
328            add_celldata_to_vtk(&mut vtk, vtkio::model::Attribute::DataArray(data_array));
329        };
330
331        sc.iter().for_each(|(r, n)| {
332            let ve: Vec<f64> = r.values.iter().map(|i| i.value).collect();
333            add_cell(n, ve)
334        });
335
336        let volumes_data = self.model.get_real_volume();
337        add_cell("real_volume", volumes_data);
338
339        let mut vtk_bytes = Vec::<u8>::new();
340        vtk.write_xml(&mut vtk_bytes).unwrap();
341        std::fs::write(path, vtk_bytes).unwrap();
342
343        Ok(())
344    }
345}
346
347#[cfg(test)]
348mod test {
349    use super::*;
350
351    ///An unreadable geometry used to be reported as "Arc error" by the io thread
352    #[test]
353    fn test_init_reports_unreadable_geometry() {
354        let error = match CMHandle::init(
355            [1, 1, 1],
356            "/nonexistent",
357            "geometry.geo",
358            grid::MeshType::Cylindrical,
359        ) {
360            Ok(_) => panic!("a missing geometry must not build a handle"),
361            Err(error) => error,
362        };
363
364        assert!(matches!(error, CoreError::IO(_)), "{}", error);
365    }
366}
cmtool_core/lib.rs

cmtool_core/
lib.rs
