honeycomb_kernels/remeshing/
capture.rs

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

honeycomb_kernels/remeshing/
capture.rs
