visualization/
result_probe.rs

1//! Read-only picking against the same deformed boundary triangles as contours.
2use crate::ContourSettings;
3use bevy::prelude::*;
4use fem_core::{Aabb, Bvh, ElementId, FemMesh, NodeId, ResultField, StepResult};
5
6/// Shared by contour rendering and probing; scale only affects the location,
7/// never the reported result value.
8pub(crate) fn deformed_position(
9    position: Vec3,
10    index: usize,
11    displacement: Option<&ResultField>,
12    scale: f32,
13) -> Vec3 {
14    if let Some(ResultField::NodeVector { values, .. }) = displacement {
15        if let Some(value) = values.get(index) {
16            return position + *value * scale;
17        }
18    }
19    position
20}
21
22struct Triangle {
23    points: [Vec3; 3],
24    nodes: [(NodeId, usize); 3],
25    element: Option<(ElementId, usize)>,
26}
27
28pub struct ProbeSurface {
29    triangles: Vec<Triangle>,
30    bvh: Bvh,
31}
32
33#[derive(Debug, Clone)]
34pub struct ProbeHit {
35    pub distance: f32,
36    /// Sample location: selected vertex for a nodal field, hit point for an element field.
37    pub point: Vec3,
38    pub element: Option<ElementId>,
39    /// Some only for nodal fields. The value is sampled at this triangle vertex,
40    /// not interpolated at the cursor and not converted to an element result.
41    pub node: Option<NodeId>,
42    pub value: Option<f32>,
43}
44
45impl ProbeSurface {
46    pub fn build(mesh: &FemMesh, step: &StepResult, settings: &ContourSettings) -> Self {
47        let mut triangles = Vec::new();
48        let Some(field) = step.field_by_name(&settings.field_name) else {
49            return Self::empty();
50        };
51        let valid = match field {
52            ResultField::NodeScalar { values, .. } => values.len() == mesh.nodes.len(),
53            ResultField::NodeVector { values, .. } => values.len() == mesh.nodes.len(),
54            ResultField::ElementScalar { values, .. } => values.len() == mesh.elements.len(),
55        };
56        if !valid {
57            return Self::empty();
58        }
59        let nodes: std::collections::HashMap<_, _> = mesh
60            .nodes
61            .iter()
62            .enumerate()
63            .map(|(i, n)| (n.id, i))
64            .collect();
65        let elements: std::collections::HashMap<_, _> = mesh
66            .elements
67            .iter()
68            .enumerate()
69            .map(|(i, e)| (e.id, i))
70            .collect();
71        let displacement = if settings.show_deformation {
72            step.field_by_name(&settings.displacement_field)
73        } else {
74            None
75        };
76        for face in mesh.cached_boundary_faces() {
77            let element = face
78                .element
79                .and_then(|id| elements.get(&id).map(|&i| (id, i)));
80            if matches!(field, ResultField::ElementScalar { .. }) && element.is_none() {
81                continue;
82            }
83            let Some(indices): Option<Vec<_>> =
84                face.nodes.iter().map(|id| nodes.get(id).copied()).collect()
85            else {
86                continue;
87            };
88            let points: Vec<_> = indices
89                .iter()
90                .map(|&i| {
91                    deformed_position(
92                        mesh.nodes[i].position,
93                        i,
94                        displacement,
95                        settings.deformation_scale,
96                    )
97                })
98                .collect();
99            if points.len() < 3 || points.iter().any(|p| !p.is_finite()) {
100                continue;
101            }
102            // Exactly the fan triangulation used by build_contour_surface_mesh.
103            for i in 1..points.len() - 1 {
104                let corners = [0, i, i + 1];
105                triangles.push(Triangle {
106                    points: corners.map(|j| points[j]),
107                    nodes: corners.map(|j| (mesh.nodes[indices[j]].id, indices[j])),
108                    element,
109                });
110            }
111        }
112        let bounds: Vec<_> = triangles
113            .iter()
114            .map(|t| {
115                let b = Aabb::from_points(&t.points).unwrap();
116                // Give flat boxes thickness, avoiding 0 * infinity in slab tests
117                // for axis-aligned camera rays exactly on a boundary.
118                b.expanded(b.extent().max_element().max(1.0) * 1e-6)
119            })
120            .collect();
121        Self {
122            bvh: Bvh::build(&bounds),
123            triangles,
124        }
125    }
126
127    fn empty() -> Self {
128        Self {
129            triangles: vec![],
130            bvh: Bvh::build(&[]),
131        }
132    }
133
134    pub fn pick(&self, origin: Vec3, direction: Vec3, field: &ResultField) -> Option<ProbeHit> {
135        if !origin.is_finite() || !direction.is_finite() {
136            return None;
137        }
138        let mut best: Option<ProbeHit> = None;
139        for index in self.bvh.query_ray(origin, direction) {
140            let triangle = &self.triangles[index as usize];
141            let Some(distance) = intersect(origin, direction, triangle.points) else {
142                continue;
143            };
144            if best.as_ref().is_some_and(|h| h.distance <= distance) {
145                continue;
146            }
147            let point = origin + direction * distance;
148            let corner = (0..3)
149                .min_by(|&a, &b| {
150                    triangle.points[a]
151                        .distance_squared(point)
152                        .total_cmp(&triangle.points[b].distance_squared(point))
153                })
154                .unwrap();
155            let (id, node_index) = triangle.nodes[corner];
156            let (node, value) = match field {
157                ResultField::NodeScalar { values, .. } => {
158                    (Some(id), values.get(node_index).copied())
159                }
160                ResultField::NodeVector { values, .. } => {
161                    (Some(id), values.get(node_index).map(|v| v.length()))
162                }
163                ResultField::ElementScalar { values, .. } => (
164                    None,
165                    triangle.element.and_then(|(_, i)| values.get(i).copied()),
166                ),
167            };
168            best = Some(ProbeHit {
169                distance,
170                point: if node.is_some() {
171                    triangle.points[corner]
172                } else {
173                    point
174                },
175                node,
176                value: value.filter(|v| v.is_finite()),
177                element: triangle.element.map(|(id, _)| id),
178            });
179        }
180        best
181    }
182}
183
184fn intersect(origin: Vec3, direction: Vec3, points: [Vec3; 3]) -> Option<f32> {
185    let [a, b, c] = points.map(|p| p.as_dvec3());
186    let e1 = b - a;
187    let e2 = c - a;
188    let p = direction.as_dvec3().cross(e2);
189    let determinant = e1.dot(p);
190    if determinant.abs() <= f64::EPSILON * e1.length() * e2.length() * 16.0 {
191        return None;
192    }
193    let t = origin.as_dvec3() - a;
194    let u = t.dot(p) / determinant;
195    let q = t.cross(e1);
196    let v = direction.as_dvec3().dot(q) / determinant;
197    if u < -1e-8 || v < -1e-8 || u + v > 1.0 + 1e-8 {
198        return None;
199    }
200    let distance = (e2.dot(q) / determinant) as f32;
201    (distance.is_finite() && distance > 0.0).then_some(distance)
202}
203
204#[cfg(test)]
205mod tests {
206    use super::*;
207    fn settings() -> ContourSettings {
208        ContourSettings {
209            mesh_index: 0,
210            step_index: 0,
211            field_name: "E".into(),
212            show_deformation: true,
213            displacement_field: "Displacement".into(),
214            deformation_scale: 10.0,
215        }
216    }
217    #[test]
218    fn probe_tracks_deformation_and_preserves_element_value_and_id() {
219        let mesh = FemMesh::demo_hex8();
220        let step = StepResult {
221            fields: vec![
222                ResultField::ElementScalar {
223                    name: "E".into(),
224                    values: vec![42.0],
225                    min: 42.0,
226                    max: 42.0,
227                },
228                ResultField::NodeVector {
229                    name: "Displacement".into(),
230                    values: vec![Vec3::X; 8],
231                    min_mag: 1.0,
232                    max_mag: 1.0,
233                },
234            ],
235            ..default()
236        };
237        let center = mesh.nodes.iter().map(|n| n.position).sum::<Vec3>() / mesh.nodes.len() as f32;
238        let index = ProbeSurface::build(&mesh, &step, &settings());
239        let field = &step.fields[0];
240        assert!(
241            index
242                .pick(center + Vec3::Z * 100.0, Vec3::NEG_Z, field)
243                .is_none()
244        );
245        let hit = index
246            .pick(
247                center + Vec3::X * 10.0 + Vec3::Z * 100.0,
248                Vec3::NEG_Z,
249                field,
250            )
251            .unwrap();
252        assert_eq!(hit.element, Some(mesh.elements[0].id));
253        assert_eq!(hit.node, None);
254        assert_eq!(hit.value, Some(42.0));
255        let mut off = settings();
256        off.show_deformation = false;
257        assert!(
258            ProbeSurface::build(&mesh, &step, &off)
259                .pick(center + Vec3::Z * 100.0, Vec3::NEG_Z, field)
260                .is_some()
261        );
262    }
263    #[test]
264    fn nodal_values_are_not_interpolated_or_zero_filled() {
265        let mesh = FemMesh::demo_hex8();
266        let mut settings = settings();
267        settings.show_deformation = false;
268        settings.field_name = "N".into();
269        let mut field = ResultField::NodeScalar {
270            name: "N".into(),
271            values: (0..8).map(|i| i as f32 + 10.0).collect(),
272            min: 10.0,
273            max: 17.0,
274        };
275        let step = StepResult {
276            fields: vec![field.clone()],
277            ..default()
278        };
279        let index = ProbeSurface::build(&mesh, &step, &settings);
280        let center = mesh.nodes.iter().map(|n| n.position).sum::<Vec3>() / 8.0;
281        let origin = center + Vec3::Z * 100.0;
282        let hit = index.pick(origin, Vec3::NEG_Z, &field).unwrap();
283        let i = mesh
284            .nodes
285            .iter()
286            .position(|n| Some(n.id) == hit.node)
287            .unwrap();
288        assert_eq!(hit.value, Some(i as f32 + 10.0));
289        if let ResultField::NodeScalar { values, .. } = &mut field {
290            values.fill(f32::NAN);
291        }
292        assert_eq!(index.pick(origin, Vec3::NEG_Z, &field).unwrap().value, None);
293    }
294
295    #[test]
296    fn probe_triangles_match_rendered_positions_and_pick_front_surface() {
297        use bevy::mesh::{Indices, VertexAttributeValues};
298        let mesh = FemMesh::demo_hex8();
299        let mut settings = settings();
300        settings.field_name = "Displacement".into();
301        let step = StepResult {
302            fields: vec![ResultField::NodeVector {
303                name: "Displacement".into(),
304                values: vec![Vec3::new(0.3, 0.4, 0.0); mesh.nodes.len()],
305                min_mag: 0.5,
306                max_mag: 0.5,
307            }],
308            ..default()
309        };
310        let rendered =
311            crate::demo_mesh::build_contour_surface_mesh(&mesh, &step, &settings, None).unwrap();
312        let Some(VertexAttributeValues::Float32x3(vertices)) =
313            rendered.attribute(Mesh::ATTRIBUTE_POSITION)
314        else {
315            panic!("expected surface positions");
316        };
317        let indices: Vec<usize> = match rendered.indices() {
318            Some(Indices::U16(v)) => v.iter().map(|&i| i as usize).collect(),
319            Some(Indices::U32(v)) => v.iter().map(|&i| i as usize).collect(),
320            None => (0..vertices.len()).collect(),
321        };
322        let surface = ProbeSurface::build(&mesh, &step, &settings);
323        assert_eq!(surface.triangles.len() * 3, indices.len());
324        for (triangle, indices) in surface.triangles.iter().zip(indices.chunks_exact(3)) {
325            for (p, &i) in triangle.points.iter().zip(indices) {
326                assert_eq!(*p, Vec3::from_array(vertices[i]));
327            }
328        }
329        let center = mesh.nodes.iter().map(|n| n.position).sum::<Vec3>() / mesh.nodes.len() as f32;
330        let hit = surface
331            .pick(
332                center + Vec3::new(3.0, 4.0, 100.0),
333                Vec3::NEG_Z,
334                &step.fields[0],
335            )
336            .unwrap();
337        let front_z = mesh
338            .nodes
339            .iter()
340            .map(|n| n.position.z)
341            .fold(f32::NEG_INFINITY, f32::max);
342        assert_eq!(hit.point.z, front_z);
343        assert!((hit.distance - (center.z + 100.0 - front_z)).abs() < 1e-5);
344        assert_eq!(hit.value, Some(0.5)); // magnitude, not multiplied by deformation scale
345    }
346}
visualization/result_probe.rs

visualization/
result_probe.rs
