honeycomb_kernels/remeshing/
capture.rs

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

honeycomb_kernels/remeshing/
capture.rs
