honeycomb_kernels/remeshing/
capture.rs

1use std::collections::{HashSet, VecDeque};
2
3use honeycomb_core::{
4    cmap::{CMap2, CMapBuilder, DartIdType, GridDescriptor, NULL_DART_ID, OrbitPolicy},
5    geometry::CoordsFloat,
6};
7use vtkio::Vtk;
8
9use crate::grisubal::{
10    Clip, GrisubalError,
11    model::{Boundary, Geometry2},
12    routines::{
13        clip_left, clip_right, compute_intersection_ids, compute_overlapping_grid,
14        detect_orientation_issue, generate_edge_data, generate_intersection_data,
15        group_intersections_per_edge, insert_edges_in_map, insert_intersections,
16    },
17};
18use crate::utils::{CurveIdType, EdgeAnchor, FaceAnchor, VertexAnchor};
19
20#[allow(clippy::missing_errors_doc)]
21/// Capture the geometry specified as input using the `grisubal` algorithm.
22///
23/// This function follows roughly the same step as the `grisubal` algorithm, but differs in two
24/// instances:
25/// - The overlapping grid avoids landing Points of Interest on all edges of the grid (instead of
26///   only corners), meaning we do not remove any from the original geometry.
27/// - Points of Interest are used to associate `VertexAnchor::Node` values to some vertices of the mesh.
28///
29/// A complete classification of the resulting mesh can be obtained using the [`classify_capture`] function.
30///
31/// # Return / Errors
32///
33/// This function returns a `Result` taking the following values:
34/// - `Ok(CMap2)` -- Algorithm ran successfully.
35/// - `Err(GrisubalError)` -- Algorithm encountered an issue. See [`GrisubalError`] for all
36///   possible errors.
37///
38/// # Panics
39///
40/// This function may panic if the specified VTK file cannot be opened.
41#[allow(clippy::needless_pass_by_value)]
42pub fn capture_geometry<T: CoordsFloat>(
43    file_path: impl AsRef<std::path::Path>,
44    grid_cell_sizes: [T; 2],
45    clip: Clip,
46) -> Result<CMap2<T>, GrisubalError> {
47    // --- IMPORT VTK INPUT
48    let geometry_vtk = match Vtk::import(file_path) {
49        Ok(vtk) => vtk,
50        Err(e) => panic!("E: could not open specified vtk file - {e}"),
51    };
52
53    // --- BUILD OUR MODEL FROM THE VTK IMPORT
54    let geometry = Geometry2::try_from(geometry_vtk)?;
55
56    // --- FIRST DETECTION OF ORIENTATION ISSUES
57    detect_orientation_issue(&geometry)?;
58
59    // --- FIND AN OVERLAPPING GRID
60    let ([nx, ny], origin) = compute_overlapping_grid(&geometry, grid_cell_sizes, true)?;
61    let [cx, cy] = grid_cell_sizes;
62    let ogrid = GridDescriptor::default()
63        .n_cells([nx, ny])
64        .len_per_cell([cx, cy])
65        .origin([origin.0, origin.1]);
66
67    // --- BUILD THE GRID
68    let mut cmap = CMapBuilder::from_grid_descriptor(ogrid)
69        .add_attribute::<Boundary>() // will be used for clipping
70        .add_attribute::<VertexAnchor>()
71        .add_attribute::<EdgeAnchor>()
72        .add_attribute::<FaceAnchor>()
73        .build()
74        .expect("E: unreachable"); // unreachable because grid dims are valid
75
76    // process the geometry
77
78    // --- STEP 1 & 2
79    // (1)
80    let (new_segments, intersection_metadata) =
81        generate_intersection_data(&cmap, &geometry, [nx, ny], [cx, cy], origin);
82    // (2)
83    let n_intersec = intersection_metadata.len();
84    let (edge_intersec, dart_slices) =
85        group_intersections_per_edge(&mut cmap, intersection_metadata);
86    let intersection_darts = compute_intersection_ids(n_intersec, &edge_intersec, &dart_slices);
87
88    // --- STEP 3
89    insert_intersections(&cmap, &edge_intersec, &dart_slices);
90
91    // --- STEP 4
92    let edges = generate_edge_data(&cmap, &geometry, &new_segments, &intersection_darts);
93
94    // --- STEP 5
95    insert_edges_in_map(&mut cmap, &edges);
96
97    // --- CLIP
98    match clip {
99        Clip::Left => clip_left(&mut cmap)?,
100        Clip::Right => clip_right(&mut cmap)?,
101        Clip::None => {}
102    }
103
104    // CLEANUP
105    cmap.remove_attribute_storage::<Boundary>();
106
107    Ok(cmap)
108}
109
110/// Error-modeling enum for classification issues.
111#[derive(thiserror::Error, Debug, PartialEq, Eq)]
112pub enum ClassificationError {
113    /// One of the attribute used to classify mesh entities isn't included in the map.
114    #[error("missing attribute for classification: {0}")]
115    MissingAttribute(&'static str),
116    /// There is an unsupported configuration in the geometry.
117    #[error("geometry contains an ambiguous or unsupported configuration: {0}")]
118    UnsupportedGeometry(&'static str),
119}
120
121/// Classify all entities of a mesh on a default geometry.
122///
123/// This function classifies all i-cells of a mesh using a basic traversal algorithm. It expects
124/// the map passed as argument to already have a set of vertices linked to nodes.
125///
126/// The algorithm uses a first oriented traversal starting from anchored vertices to classifies
127/// boundaries as curves. It also checks for curves that are not linked to anchored vertices (e.g.
128/// a hole inside of the geometry).
129///
130/// Once curves are classified, the remaining unclassified cells are associated to surfaces, which
131/// are identified using a BFS-like algorithm stopping on already marked edges, to correctly
132/// discriminate surfaces.
133///
134/// # Errors
135///
136/// This function may fail and return an error if:
137/// - one of the attribute used to classify entities is missing (e.g. `VertexAnchor`),
138/// - The structure of the mesh is incorrect / unsupported (e.g. an open geometry).
139///
140/// # Panics
141///
142/// In `debug` mode, we use assertions to check every single i-cell of the mesh has been
143/// classified. If that is not the case, the function fill panic.
144#[allow(clippy::cast_possible_truncation, clippy::too_many_lines)]
145pub fn classify_capture<T: CoordsFloat>(cmap: &CMap2<T>) -> Result<(), ClassificationError> {
146    if !cmap.contains_attribute::<VertexAnchor>() {
147        Err(ClassificationError::MissingAttribute(
148            std::any::type_name::<VertexAnchor>(),
149        ))?;
150    }
151    if !cmap.contains_attribute::<EdgeAnchor>() {
152        Err(ClassificationError::MissingAttribute(
153            std::any::type_name::<EdgeAnchor>(),
154        ))?;
155    }
156    if !cmap.contains_attribute::<FaceAnchor>() {
157        Err(ClassificationError::MissingAttribute(
158            std::any::type_name::<FaceAnchor>(),
159        ))?;
160    }
161
162    let mut curve_id = 0;
163
164    // classify boundaries
165    for (i, dart) in cmap
166        .iter_vertices()
167        .filter_map(|v| {
168            cmap.force_read_attribute::<VertexAnchor>(v).and_then(|_| {
169                cmap.orbit(OrbitPolicy::Vertex, v)
170                    .find(|d| cmap.beta::<2>(*d) == NULL_DART_ID)
171            })
172        })
173        .enumerate()
174    {
175        mark_curve(cmap, dart, i as CurveIdType)?;
176        curve_id = curve_id.max(i);
177    }
178    // check for boundaries that are not reachable from anchored vertices (e.g. a hole with no PoI)
179    while let Some(dart) = (1..cmap.n_darts() as DartIdType)
180        .filter_map(|d| {
181            // check only used darts on the boundary
182            let used = !cmap.is_unused(d);
183            if used {
184                cmap.orbit(OrbitPolicy::Vertex, d)
185                    .find(|dd| cmap.beta::<2>(*dd) == NULL_DART_ID)
186            } else {
187                None
188            }
189        })
190        .find(|d| {
191            cmap.force_read_attribute::<EdgeAnchor>(cmap.edge_id(*d))
192                .is_none()
193        })
194    {
195        curve_id += 1; // use a new curve id
196        cmap.force_write_attribute(
197            cmap.vertex_id(dart),
198            VertexAnchor::Curve(curve_id as CurveIdType),
199        );
200        mark_curve(cmap, dart, curve_id as CurveIdType)?;
201    }
202
203    // classify inner entities using a coloring-like algorithm
204    let mut surface_id = 0;
205    let mut queue = VecDeque::new();
206    let mut marked = HashSet::new();
207    marked.insert(0);
208    cmap.iter_faces()
209        // does this filter item updated in the for_each block?
210        // .filter(|f| cmap.force_read_attribute::<FaceAnchor>(*f).is_some())
211        .for_each(|f| {
212            if cmap.force_read_attribute::<FaceAnchor>(f).is_none() {
213                queue.clear();
214                queue.push_front(f);
215                while let Some(crt) = queue.pop_front() {
216                    // if the filter works correctly, this if isn't useful
217                    cmap.force_write_attribute(crt, FaceAnchor::Surface(surface_id));
218                    cmap.orbit(OrbitPolicy::Face, crt as DartIdType)
219                        .filter(|d| {
220                            cmap.force_read_attribute::<EdgeAnchor>(cmap.edge_id(*d))
221                                .is_none()
222                        })
223                        .for_each(|d| {
224                            cmap.force_write_attribute(
225                                cmap.edge_id(d),
226                                EdgeAnchor::Surface(surface_id),
227                            );
228                            if cmap
229                                .force_read_attribute::<VertexAnchor>(cmap.vertex_id(d))
230                                .is_none()
231                            {
232                                cmap.force_write_attribute(
233                                    cmap.vertex_id(d),
234                                    VertexAnchor::Surface(surface_id),
235                                );
236                            }
237                            let neighbor_face = cmap.face_id(cmap.beta::<2>(d));
238                            if marked.insert(neighbor_face) {
239                                queue.push_back(neighbor_face);
240                            }
241                        });
242                }
243                surface_id += 1;
244            }
245        });
246
247    // in debug mode, ensure all entities are classified
248    debug_assert!(
249        cmap.iter_vertices()
250            .all(|v| cmap.force_read_attribute::<VertexAnchor>(v).is_some()),
251        "E: Not all vertices are classified",
252    );
253    debug_assert!(
254        cmap.iter_edges()
255            .all(|e| cmap.force_read_attribute::<EdgeAnchor>(e).is_some()),
256        "E: Not all edges are classified",
257    );
258    debug_assert!(
259        cmap.iter_faces()
260            .all(|f| cmap.force_read_attribute::<FaceAnchor>(f).is_some()),
261        "E: Not all faces are classified",
262    );
263
264    Ok(())
265}
266
267/// Traverse and anchor a boundary starting from the specified dart.
268///
269/// All entities making up the boundary (vertices, edges) are anchored to the `curve_id` curve.
270/// The only exception is the vertex associated to the starting dart.
271#[allow(clippy::cast_possible_truncation)]
272fn mark_curve<T: CoordsFloat>(
273    cmap: &CMap2<T>,
274    start: DartIdType,
275    curve_id: CurveIdType,
276) -> Result<(), ClassificationError> {
277    // only write attribute on the edge for the first one
278    // since we start from an anchored vertex
279    cmap.force_write_attribute::<EdgeAnchor>(cmap.edge_id(start), EdgeAnchor::Curve(curve_id));
280    let mut next = cmap.beta::<1>(start);
281    while cmap
282        .force_read_attribute::<VertexAnchor>(cmap.vertex_id(next))
283        .is_none()
284    {
285        if let Some(crt) = cmap
286            .orbit(OrbitPolicy::Vertex, next)
287            .find(|d| cmap.beta::<2>(*d) == NULL_DART_ID)
288        {
289            cmap.force_write_attribute(cmap.vertex_id(crt), VertexAnchor::Curve(curve_id));
290            cmap.force_write_attribute(cmap.edge_id(crt), EdgeAnchor::Curve(curve_id));
291            next = cmap.beta::<1>(crt);
292        } else {
293            // this should be unreachable as long as the geometry is closed and the node is on the boundary
294            Err(ClassificationError::UnsupportedGeometry(
295                "open geometry or node outside of boundary",
296            ))?;
297        }
298    }
299    Ok(())
300}
honeycomb_kernels/remeshing/capture.rs

honeycomb_kernels/remeshing/
capture.rs
