diff --git a/.cargo/config.toml b/.cargo/config.toml new file mode 100644 index 0000000..e11d56a --- /dev/null +++ b/.cargo/config.toml @@ -0,0 +1,2 @@ +[alias] +xtask = "run --release --package xtask --" diff --git a/.gitignore b/.gitignore index d67dbc0..76a071b 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ build/ffi/ *.o *.a target/ +assets/old/ Cargo.lock diff --git a/Cargo.toml b/Cargo.toml index 098cb69..c0e6725 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,7 @@ [workspace] resolver = "2" -members = ["crates/*"] +members = ["crates/*", "xtask"] +default-members = ["crates/femm-app"] [workspace.package] edition = "2024" diff --git a/assets/belasolv.svg b/assets/belasolv.svg new file mode 100644 index 0000000..5a8ad0c --- /dev/null +++ b/assets/belasolv.svg @@ -0,0 +1,9 @@ + + + + + + + + + \ No newline at end of file diff --git a/assets/femm.svg b/assets/femm.svg new file mode 100644 index 0000000..3f24543 --- /dev/null +++ b/assets/femm.svg @@ -0,0 +1,5 @@ + + + + + \ No newline at end of file diff --git a/assets/idr_cviewtype.svg b/assets/idr_cviewtype.svg new file mode 100644 index 0000000..bf23ee9 --- /dev/null +++ b/assets/idr_cviewtype.svg @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/assets/idr_femmplottype.svg b/assets/idr_femmplottype.svg new file mode 100644 index 0000000..c86a47f --- /dev/null +++ b/assets/idr_femmplottype.svg @@ -0,0 +1,7 @@ + + + + + + + \ No newline at end of file diff --git a/crates/femm-app/src/doc_canvas.rs b/crates/femm-app/src/doc_canvas.rs index 7b8065f..420b7ec 100644 --- a/crates/femm-app/src/doc_canvas.rs +++ b/crates/femm-app/src/doc_canvas.rs @@ -18,6 +18,8 @@ const BG: Color = Color::WHITE; const GEOM: Color = Color::BLACK; const LABEL_COLOR: Color = Color::from_rgb(0.25, 0.45, 0.85); const PENDING_COLOR: Color = Color::from_rgb(0.85, 0.30, 0.30); +const SELECT_COLOR: Color = Color::from_rgb(0.95, 0.20, 0.20); +const SELECT_STROKE: f32 = 2.4; /// active editing mode on the canvas. #[derive(Debug, Default, Clone, Copy, PartialEq, Eq)] @@ -36,6 +38,10 @@ pub enum CanvasMessage { Click { world: (f64, f64), tool: Tool }, /// two-point segment request from the canvas. SegmentBetween { from: (f64, f64), to: (f64, f64) }, + /// right-click on Select mode, toggling the closest entity to the given world point. + TogglePickAt { world: (f64, f64) }, + /// Delete key in Select mode, removing every selected entity from the doc. + DeleteSelected, } /// pan offset, zoom factor, and click-gesture bookkeeping for the canvas. @@ -84,6 +90,16 @@ impl<'a> canvas::Program for DocCanvas<'a> { return Some(Action::capture()); } } + Event::Mouse(mouse::Event::ButtonReleased(mouse::Button::Right)) => { + if self.tool == Tool::Select { + if let Some(now) = cursor.position_in(bounds) { + let view = ViewTransform::fit(self.doc, bounds, state); + let world = view.inverse_map(now); + return Some(Action::publish(CanvasMessage::TogglePickAt { world }) + .and_capture()); + } + } + } Event::Mouse(mouse::Event::ButtonReleased(mouse::Button::Left)) => { let press = state.press_origin.take(); state.drag_origin = None; @@ -162,11 +178,23 @@ impl<'a> canvas::Program for DocCanvas<'a> { return Some(Action::request_redraw()); } } - if matches!(key, iced::keyboard::Key::Named(iced::keyboard::key::Named::Escape)) { + if matches!( + key, + iced::keyboard::Key::Named(iced::keyboard::key::Named::Escape), + ) { if state.pending_segment_start.take().is_some() { return Some(Action::request_redraw()); } } + if matches!( + key, + iced::keyboard::Key::Named(iced::keyboard::key::Named::Delete) + | iced::keyboard::Key::Named(iced::keyboard::key::Named::Backspace), + ) { + if self.tool == Tool::Select { + return Some(Action::publish(CanvasMessage::DeleteSelected)); + } + } } _ => {} } @@ -192,8 +220,13 @@ impl<'a> canvas::Program for DocCanvas<'a> { { let a = view.map(p0.x, p0.y); let b = view.map(p1.x, p1.y); + let (color, width) = if s.selected { + (SELECT_COLOR, SELECT_STROKE) + } else { + (GEOM, STROKE_WIDTH) + }; frame.stroke(&Path::line(a, b), - Stroke::default().with_width(STROKE_WIDTH).with_color(GEOM)); + Stroke::default().with_width(width).with_color(color)); } } @@ -201,6 +234,11 @@ impl<'a> canvas::Program for DocCanvas<'a> { if let (Some(p0), Some(p1)) = (self.doc.nodes.get(a.n0 as usize), self.doc.nodes.get(a.n1 as usize)) { + let (color, width) = if a.selected { + (SELECT_COLOR, SELECT_STROKE) + } else { + (GEOM, STROKE_WIDTH) + }; if let Some((center, radius, start_angle, end_angle)) = arc_geometry(p0.x, p0.y, p1.x, p1.y, a.arc_length, a.normal_direction) { @@ -214,19 +252,24 @@ impl<'a> canvas::Program for DocCanvas<'a> { end_angle, }); frame.stroke(&b.build(), - Stroke::default().with_width(STROKE_WIDTH).with_color(GEOM)); + Stroke::default().with_width(width).with_color(color)); } else { let a_px = view.map(p0.x, p0.y); let b_px = view.map(p1.x, p1.y); frame.stroke(&Path::line(a_px, b_px), - Stroke::default().with_width(STROKE_WIDTH).with_color(GEOM)); + Stroke::default().with_width(width).with_color(color)); } } } for n in &self.doc.nodes { let p = view.map(n.x, n.y); - frame.fill(&Path::circle(p, NODE_RADIUS), GEOM); + let (color, r) = if n.selected { + (SELECT_COLOR, NODE_RADIUS + 2.0) + } else { + (GEOM, NODE_RADIUS) + }; + frame.fill(&Path::circle(p, r), color); } if let Some(start_world) = state.pending_segment_start { @@ -245,6 +288,11 @@ impl<'a> canvas::Program for DocCanvas<'a> { for label in &self.doc.block_labels { let p = view.map(label.x, label.y); + let (color, width) = if label.selected { + (SELECT_COLOR, SELECT_STROKE) + } else { + (LABEL_COLOR, STROKE_WIDTH) + }; let cross = Path::new(|b| { b.move_to(Point::new(p.x - LABEL_TICK_PX, p.y)); b.line_to(Point::new(p.x + LABEL_TICK_PX, p.y)); @@ -252,11 +300,11 @@ impl<'a> canvas::Program for DocCanvas<'a> { b.line_to(Point::new(p.x, p.y + LABEL_TICK_PX)); }); frame.stroke(&cross, - Stroke::default().with_width(STROKE_WIDTH).with_color(LABEL_COLOR)); + Stroke::default().with_width(width).with_color(color)); frame.fill_text(Text { content: label.block_type.clone(), position: Point::new(p.x + LABEL_TICK_PX + 4.0, p.y - 8.0), - color: LABEL_COLOR, + color, size: 12.0.into(), ..Text::default() }); diff --git a/crates/femm-app/src/main.rs b/crates/femm-app/src/main.rs index a5179c3..07cff7e 100644 --- a/crates/femm-app/src/main.rs +++ b/crates/femm-app/src/main.rs @@ -94,6 +94,22 @@ impl App { self.status = format!("rejected segment {n0} -> {n1}"); } } + Message::Canvas(CanvasMessage::TogglePickAt { world }) => { + toggle_closest(&mut self.doc, world.0, world.1); + self.status = String::from("toggled selection at right-click"); + } + Message::Canvas(CanvasMessage::DeleteSelected) => { + let n = self.doc.delete_selected_nodes(); + let s = self.doc.delete_selected_segments(); + let a = self.doc.delete_selected_arcs(); + let b = self.doc.delete_selected_block_labels(); + let total = n + s + a + b; + self.status = if total == 0 { + String::from("nothing selected") + } else { + format!("deleted: {n} nodes, {s} segments, {a} arcs, {b} labels") + }; + } } Task::none() } @@ -134,6 +150,54 @@ impl App { } } +/// flips the selected flag on whichever entity sits nearest to (x, y) in doc world coords. +fn toggle_closest(doc: &mut FemmDoc, x: f64, y: f64) { + use femm_doc_mag::geom_math::{ + shortest_distance_from_arc, shortest_distance_from_segment, + }; + + #[derive(Clone, Copy)] + enum Kind { Node, Segment, Arc, Label } + let mut best: Option<(Kind, usize, f64)> = None; + let mut consider = |kind: Kind, idx: usize, d: f64| { + match best { + None => best = Some((kind, idx, d)), + Some((_, _, bd)) if d < bd => best = Some((kind, idx, d)), + _ => {} + } + }; + + for (i, n) in doc.nodes.iter().enumerate() { + consider(Kind::Node, i, (n.x - x).hypot(n.y - y)); + } + for (i, s) in doc.segments.iter().enumerate() { + if let (Some(p0), Some(p1)) = + (doc.nodes.get(s.n0 as usize), doc.nodes.get(s.n1 as usize)) + { + consider(Kind::Segment, i, shortest_distance_from_segment((x, y), (p0.x, p0.y), (p1.x, p1.y))); + } + } + for (i, a) in doc.arcs.iter().enumerate() { + if let (Some(p0), Some(p1)) = + (doc.nodes.get(a.n0 as usize), doc.nodes.get(a.n1 as usize)) + { + consider(Kind::Arc, i, shortest_distance_from_arc((x, y), (p0.x, p0.y), (p1.x, p1.y), a.arc_length)); + } + } + for (i, b) in doc.block_labels.iter().enumerate() { + consider(Kind::Label, i, (b.x - x).hypot(b.y - y)); + } + + if let Some((kind, idx, _)) = best { + match kind { + Kind::Node => doc.nodes[idx].selected ^= true, + Kind::Segment => doc.segments[idx].selected ^= true, + Kind::Arc => doc.arcs[idx].selected ^= true, + Kind::Label => doc.block_labels[idx].selected ^= true, + } + } +} + fn tool_button(label: &str, this_tool: Tool, active: Tool) -> Element<'_, Message> { let btn = button(text(label).size(13)).on_press(Message::SelectTool(this_tool)); if this_tool == active { diff --git a/crates/femm-doc-curr/src/edit.rs b/crates/femm-doc-curr/src/edit.rs new file mode 100644 index 0000000..6d773d6 --- /dev/null +++ b/crates/femm-doc-curr/src/edit.rs @@ -0,0 +1,608 @@ +//! geometry editing primitives on [`FemmDoc`]: add, delete, closest-point queries. + +use crate::geom_math::{ + arc_arc_intersection, circle_from_arc, line_arc_intersection, line_line_intersection, + shortest_distance_from_arc, shortest_distance_from_segment, +}; +use crate::{ArcSegment, BlockLabel, FemmDoc, Node, Segment}; +use num_complex::Complex64; + +/// fraction of the node-bbox diagonal used as auto-tolerance for intersection-node coalescing. +const BBOX_TOLERANCE_FRAC: f64 = 1.0e-6; +/// fraction of a segment's length within which an off-endpoint node triggers a recursive split. +const ON_LINE_FRAC: f64 = 1.0e-5; + +impl FemmDoc { + /// adds a node at (x, y), returning the index of an existing node within `tol` distance when present. + pub fn add_node(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, n) in self.nodes.iter().enumerate() { + let dx = n.x - x; + let dy = n.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.nodes.push(Node { + x, y, + boundary_marker: String::new(), + in_group: 0, + selected: false, + }); + self.nodes.len() - 1 + } + + /// adds a block label at (x, y), returning the index of an existing label within `tol` when present. + pub fn add_block_label(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, b) in self.block_labels.iter().enumerate() { + let dx = b.x - x; + let dy = b.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.block_labels.push(BlockLabel { + x, y, + max_area: 0.0, + block_type: String::from(""), + in_conductor: String::from(""), + in_group: 0, + is_external: false, + is_default: false, + selected: false, + }); + self.block_labels.len() - 1 + } + + /// adds a segment between two node indices, splitting at every crossing and through any on-line node. + pub fn add_segment(&mut self, n0: i32, n1: i32) -> bool { + self.add_segment_with_marker(n0, n1, "") + } + + /// PSLG-aware variant propagating a boundary-marker name onto every resulting segment piece. + pub fn add_segment_with_marker(&mut self, n0: i32, n1: i32, marker: &str) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // duplicate either-orientation. + let (a, b) = if n0 < n1 { (n0, n1) } else { (n1, n0) }; + for s in &self.segments { + let (sa, sb) = if s.n0 < s.n1 { (s.n0, s.n1) } else { (s.n1, s.n0) }; + if sa == a && sb == b { return false; } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + if let Some(hit) = line_line_intersection(n0p, n1p, sp0, sp1) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in line_arc_intersection(n0p, n1p, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + self.segments.push(Segment { + n0, n1, + max_side_length: -1.0, + boundary_marker: marker.to_string(), + hidden: false, + in_group: 0, + selected: false, + }); + + // first non-endpoint node on the new segment's interior, if any. + let length = ((n1p.0 - n0p.0).powi(2) + (n1p.1 - n0p.1).powi(2)).sqrt(); + let dmin = length * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_segment(np, n0p, n1p); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.segments.pop(); + let a = self.add_segment_with_marker(n0, mid, marker); + let b = self.add_segment_with_marker(mid, n1, marker); + a || b + } else { + true + } + } + + /// auto-tolerance for intersection-node coalescing, derived from the node bounding box. + fn bbox_tolerance(&self) -> f64 { + nodes_bbox_tolerance(&self.nodes) + } + + /// rebuilds every list through the PSLG-aware add primitives, catching crossings missed by incremental edits. + pub fn enforce_pslg(&mut self) { + let old_nodes = std::mem::take(&mut self.nodes); + let old_segments = std::mem::take(&mut self.segments); + let old_arcs = std::mem::take(&mut self.arcs); + let old_block_labels = std::mem::take(&mut self.block_labels); + + let tol = nodes_bbox_tolerance(&old_nodes); + + // dedupes by position with metadata preserved. + for n in &old_nodes { + let mut duplicate = false; + for existing in &self.nodes { + if (existing.x - n.x).hypot(existing.y - n.y) < tol { + duplicate = true; + break; + } + } + if !duplicate { + self.nodes.push(n.clone()); + } + } + + for s in old_segments { + let p0 = (old_nodes[s.n0 as usize].x, old_nodes[s.n0 as usize].y); + let p1 = (old_nodes[s.n1 as usize].x, old_nodes[s.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_segment_with_marker(n0 as i32, n1 as i32, &s.boundary_marker); + } + + for a in old_arcs { + let p0 = (old_nodes[a.n0 as usize].x, old_nodes[a.n0 as usize].y); + let p1 = (old_nodes[a.n1 as usize].x, old_nodes[a.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_arc_segment_with_template(n0 as i32, n1 as i32, a.arc_length, &a); + } + + self.block_labels = old_block_labels; + } + + /// adds an arc between two node indices, splitting at every crossing and through any on-arc node. + pub fn add_arc_segment(&mut self, n0: i32, n1: i32, arc_length_deg: f64) -> bool { + let template = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + max_side_length: 10.0, + boundary_marker: String::new(), + hidden: false, + in_group: 0, + normal_direction: true, + selected: false, + }; + self.add_arc_segment_with_template(n0, n1, arc_length_deg, &template) + } + + /// PSLG-aware variant propagating boundary marker and side-length metadata onto every arc piece. + pub fn add_arc_segment_with_template( + &mut self, + n0: i32, + n1: i32, + arc_length_deg: f64, + template: &ArcSegment, + ) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // same directed endpoints with similar sweep counts as duplicate. + for a in &self.arcs { + if a.n0 == n0 && a.n1 == n1 && (a.arc_length - arc_length_deg).abs() < 1.0e-2 { + return false; + } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + for hit in line_arc_intersection(sp0, sp1, n0p, n1p, arc_length_deg) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in arc_arc_intersection(n0p, n1p, arc_length_deg, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + let new_arc = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + ..template.clone() + }; + self.arcs.push(new_arc); + + // first non-endpoint node on the new arc's sweep range, if any. + let (cx, cy, radius) = circle_from_arc(n0p, n1p, arc_length_deg); + let sweep_rad = arc_length_deg.to_radians(); + let arc_length_world = radius * sweep_rad.abs(); + let dmin = arc_length_world * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_arc(np, n0p, n1p, arc_length_deg); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.arcs.pop(); + let mid_pos = (self.nodes[mid as usize].x, self.nodes[mid as usize].y); + let c = Complex64::new(cx, cy); + let a0 = Complex64::new(n0p.0, n0p.1); + let a1 = Complex64::new(n1p.0, n1p.1); + let a2 = Complex64::new(mid_pos.0, mid_pos.1); + let sweep_to_mid = ((a2 - c) / (a0 - c)).arg().to_degrees(); + let sweep_from_mid = ((a1 - c) / (a2 - c)).arg().to_degrees(); + let a = self.add_arc_segment_with_template(n0, mid, sweep_to_mid, template); + let b = self.add_arc_segment_with_template(mid, n1, sweep_from_mid, template); + a || b + } else { + true + } + } + + /// removes selected nodes and rewrites segment/arc endpoint indices to drop references. + pub fn delete_selected_nodes(&mut self) -> usize { + let keep: Vec = self.nodes.iter().map(|n| !n.selected).collect(); + let mut remap: Vec = Vec::with_capacity(keep.len()); + let mut next: i32 = 0; + for &k in &keep { + remap.push(if k { let r = next; next += 1; r } else { -1 }); + } + let removed = self.nodes.len() - next as usize; + + let mut new_nodes = Vec::with_capacity(next as usize); + for (i, n) in self.nodes.drain(..).enumerate() { + if keep[i] { new_nodes.push(n); } + } + self.nodes = new_nodes; + + self.segments.retain(|s| { + let a = s.n0 as usize; + let b = s.n1 as usize; + a < keep.len() && b < keep.len() && keep[a] && keep[b] + }); + for s in &mut self.segments { + s.n0 = remap[s.n0 as usize]; + s.n1 = remap[s.n1 as usize]; + } + + self.arcs.retain(|a| { + let i = a.n0 as usize; + let j = a.n1 as usize; + i < keep.len() && j < keep.len() && keep[i] && keep[j] + }); + for a in &mut self.arcs { + a.n0 = remap[a.n0 as usize]; + a.n1 = remap[a.n1 as usize]; + } + + removed + } + + /// removes selected segments and returns the count. + pub fn delete_selected_segments(&mut self) -> usize { + let before = self.segments.len(); + self.segments.retain(|s| !s.selected); + before - self.segments.len() + } + + /// removes selected arc segments and returns the count. + pub fn delete_selected_arcs(&mut self) -> usize { + let before = self.arcs.len(); + self.arcs.retain(|a| !a.selected); + before - self.arcs.len() + } + + /// removes selected block labels and returns the count. + pub fn delete_selected_block_labels(&mut self) -> usize { + let before = self.block_labels.len(); + self.block_labels.retain(|b| !b.selected); + before - self.block_labels.len() + } + + /// returns the index of the closest node to (x, y), or None when the node list is empty. + pub fn closest_node(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, n) in self.nodes.iter().enumerate() { + let d = (n.x - x).hypot(n.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the closest block label to (x, y), or None when none exist. + pub fn closest_block_label(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, b) in self.block_labels.iter().enumerate() { + let d = (b.x - x).hypot(b.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the segment whose nearest point is closest to (x, y). + pub fn closest_segment(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, s) in self.segments.iter().enumerate() { + let (Some(p0), Some(p1)) = (self.nodes.get(s.n0 as usize), self.nodes.get(s.n1 as usize)) else { continue }; + let d = point_to_segment_distance(x, y, p0.x, p0.y, p1.x, p1.y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// clears the selection flag on every geometric entity in the doc. + pub fn clear_selection(&mut self) { + for n in &mut self.nodes { n.selected = false; } + for s in &mut self.segments { s.selected = false; } + for a in &mut self.arcs { a.selected = false; } + for b in &mut self.block_labels { b.selected = false; } + } +} + +/// auto-tolerance derived from the bounding box of a node slice. +fn nodes_bbox_tolerance(nodes: &[Node]) -> f64 { + if nodes.len() < 2 { + return 1.0e-8; + } + let mut xmin = f64::INFINITY; + let mut xmax = f64::NEG_INFINITY; + let mut ymin = f64::INFINITY; + let mut ymax = f64::NEG_INFINITY; + for n in nodes { + xmin = xmin.min(n.x); xmax = xmax.max(n.x); + ymin = ymin.min(n.y); ymax = ymax.max(n.y); + } + let dx = xmax - xmin; + let dy = ymax - ymin; + (dx * dx + dy * dy).sqrt() * BBOX_TOLERANCE_FRAC +} + +/// Euclidean distance from (px, py) to the segment between (ax, ay) and (bx, by). +fn point_to_segment_distance(px: f64, py: f64, ax: f64, ay: f64, bx: f64, by: f64) -> f64 { + let dx = bx - ax; + let dy = by - ay; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return (px - ax).hypot(py - ay); + } + let t = (((px - ax) * dx + (py - ay) * dy) / len2).clamp(0.0, 1.0); + let cx = ax + t * dx; + let cy = ay + t * dy; + (px - cx).hypot(py - cy) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn doc_with_corners() -> FemmDoc { + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0,-1.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d + } + + #[test] + fn duplicate_segment_rejected() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert!(!d.add_segment(0, 1)); + assert!(!d.add_segment(1, 0)); + assert_eq!(d.segments.len(), 1); + } + + #[test] + fn segment_passing_through_existing_node_splits() { + // three colinear nodes (-1,0), (1,0), (0,0). adding the outer-to-outer segment + // splits at the midpoint node, producing two pieces meeting there. + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert_eq!(d.segments.len(), 2); + let touches_mid = d.segments.iter().filter(|s| s.n0 == 2 || s.n1 == 2).count(); + assert_eq!(touches_mid, 2); + } + + #[test] + fn enforce_pslg_splits_first_segment_at_late_intersection() { + // incremental edit splits only the second of two crossing segments. enforce_pslg + // back-splits the first segment through the intersection node. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + assert_eq!(d.segments.len(), 3, "pre-enforce: one whole, one split"); + + d.enforce_pslg(); + + assert_eq!(d.nodes.len(), 5, "no node added or merged by enforce"); + assert_eq!(d.segments.len(), 4, "both segments split at the origin"); + let touches_origin = d.segments.iter().filter(|s| s.n0 == 4 || s.n1 == 4).count(); + assert_eq!(touches_origin, 4); + } + + #[test] + fn enforce_pslg_preserves_segment_boundary_marker() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment_with_marker(0, 1, "outer")); + + d.enforce_pslg(); + + assert_eq!(d.segments.len(), 1); + assert_eq!(d.segments[0].boundary_marker, "outer"); + } + + #[test] + fn enforce_pslg_preserves_arc_metadata() { + // arc carrying non-default marker, side length, group, and direction. + // enforce_pslg round-trips every field intact. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let template = ArcSegment { + n0: 0, + n1: 1, + arc_length: 90.0, + max_side_length: 5.0, + boundary_marker: String::from("outer"), + hidden: true, + in_group: 7, + normal_direction: false, + selected: false, + }; + assert!(d.add_arc_segment_with_template(0, 1, 90.0, &template)); + + d.enforce_pslg(); + + assert_eq!(d.arcs.len(), 1); + let a = &d.arcs[0]; + assert_eq!(a.boundary_marker, "outer"); + assert!((a.max_side_length - 5.0).abs() < 1e-12); + assert!(a.hidden); + assert_eq!(a.in_group, 7); + assert!(!a.normal_direction); + } + + #[test] + fn duplicate_arc_rejected() { + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(!d.add_arc_segment(0, 1, 90.0)); + // reversed endpoints define a distinct arc. + assert!(d.add_arc_segment(1, 0, 90.0)); + assert_eq!(d.arcs.len(), 2); + } + + #[test] + fn add_segment_crossing_unit_quarter_arc_splits_segment() { + // unit quarter arc from (1,0) to (0,1); horizontal line at y=0.5 crosses the arc once. + // expect a fresh intersection node near (0.866, 0.5) and the segment split into two. + let mut d = FemmDoc::default(); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d.add_node(-2.0, 0.5, 0.0); + d.add_node( 2.0, 0.5, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5); + let n4 = &d.nodes[4]; + assert!((n4.x - 0.75_f64.sqrt()).abs() < 1e-6); + assert!((n4.y - 0.5).abs() < 1e-9); + assert_eq!(d.segments.len(), 2); + assert_eq!(d.arcs.len(), 1); + } + + #[test] + fn add_arc_passing_through_existing_node_splits_into_two_arcs() { + // pre-existing node at (cos45, sin45) lies on the unit quarter arc from (1,0) to (0,1). + // adding the arc detects the on-arc node and emits two sub-arcs of 45 degrees each. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let mid_x = std::f64::consts::FRAC_1_SQRT_2; + d.add_node(mid_x, mid_x, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + + assert_eq!(d.arcs.len(), 2); + let touches_mid = d.arcs.iter().filter(|a| a.n0 == 2 || a.n1 == 2).count(); + assert_eq!(touches_mid, 2); + // split sweeps sum to the parent 90 degrees, within floating slop. + let total: f64 = d.arcs.iter().map(|a| a.arc_length).sum(); + assert!((total - 90.0).abs() < 1e-6, "sum of split sweeps = {total}"); + } + + #[test] + fn crossing_second_segment_splits_at_intersection() { + // two perpendicular segments through the origin. horizontal first, vertical second. + // the vertical add inserts an origin node and self-splits; the horizontal stays whole. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5, "intersection node added"); + let new_idx = 4; + assert!((d.nodes[new_idx].x).abs() < 1e-9); + assert!((d.nodes[new_idx].y).abs() < 1e-9); + + // expected pieces: 0->1 (horizontal whole), 2->4 and 4->3 (vertical halves). + assert_eq!(d.segments.len(), 3); + let touches_new = d.segments.iter().filter(|s| s.n0 == new_idx as i32 || s.n1 == new_idx as i32).count(); + assert_eq!(touches_new, 2); + } +} diff --git a/crates/femm-doc-curr/src/geom_math.rs b/crates/femm-doc-curr/src/geom_math.rs new file mode 100644 index 0000000..37eeb54 --- /dev/null +++ b/crates/femm-doc-curr/src/geom_math.rs @@ -0,0 +1,322 @@ +//! 2D intersection math and distance queries over (x, y) pairs and degree-sweep arcs. + +use num_complex::Complex64; + +/// closeness factor for common-endpoint and small-interval rejection. +const EPS_FRAC: f64 = 1.0e-8; + +/// recovers the center and radius of the circle through an arc's endpoints with the given sweep. +pub fn circle_from_arc(p0: (f64, f64), p1: (f64, f64), sweep_deg: f64) -> (f64, f64, f64) { + let a0 = Complex64::new(p0.0, p0.1); + let a1 = Complex64::new(p1.0, p1.1); + let chord = (a1 - a0).norm(); + let t = (a1 - a0) / chord; + let tta = sweep_deg.to_radians(); + let r = chord / (2.0 * (tta / 2.0).sin()); + let c = a0 + (Complex64::new(chord / 2.0, (r * r - chord * chord / 4.0).sqrt())) * t; + (c.re, c.im, r.abs()) +} + +/// shortest distance from a point to a segment. +pub fn shortest_distance_from_segment(p: (f64, f64), a: (f64, f64), b: (f64, f64)) -> f64 { + let dx = b.0 - a.0; + let dy = b.1 - a.1; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return ((p.0 - a.0).powi(2) + (p.1 - a.1).powi(2)).sqrt(); + } + let t = (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / len2).clamp(0.0, 1.0); + let cx = a.0 + t * dx; + let cy = a.1 + t * dy; + ((p.0 - cx).powi(2) + (p.1 - cy).powi(2)).sqrt() +} + +/// shortest distance from a point to a sweep-defined arc. +pub fn shortest_distance_from_arc( + p: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> f64 { + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let pc = Complex64::new(p.0 - cx, p.1 - cy); + let d = pc.norm(); + if d == 0.0 { + return r; + } + let t = pc / d; + let foot = pc - Complex64::new(r * t.re, r * t.im); + let l = foot.norm(); + let a0_dir = Complex64::new(a0.0 - cx, a0.1 - cy); + let z = ((t / a0_dir).arg() * 180.0 / std::f64::consts::PI + 360.0) % 360.0; + let sweep_abs = sweep_deg.abs(); + if z > 0.0 && z < sweep_abs { + return l; + } + let e0 = ((p.0 - a0.0).powi(2) + (p.1 - a0.1).powi(2)).sqrt(); + let e1 = ((p.0 - a1.0).powi(2) + (p.1 - a1.1).powi(2)).sqrt(); + e0.min(e1) +} + +/// intersection of the open segments p0->p1 and q0->q1, or None when none exists. +pub fn line_line_intersection( + p0: (f64, f64), + p1: (f64, f64), + q0: (f64, f64), + q1: (f64, f64), +) -> Option<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let q0c = Complex64::new(q0.0, q0.1); + let q1c = Complex64::new(q1.0, q1.1); + + // shared endpoint -> no other intersection + if (p0c - q0c).norm() < f64::EPSILON + || (p0c - q1c).norm() < f64::EPSILON + || (p1c - q0c).norm() < f64::EPSILON + || (p1c - q1c).norm() < f64::EPSILON + { + return None; + } + + let ee = ((p1c - p0c).norm()).min((q1c - q0c).norm()) * EPS_FRAC; + + let denom = p1c - p0c; + if denom.norm() < f64::EPSILON { + return None; + } + let r0 = (q0c - p0c) / denom; + let r1 = (q1c - p0c) / denom; + + if r0.re <= 0.0 && r1.re <= 0.0 { return None; } + if r0.re >= 1.0 && r1.re >= 1.0 { return None; } + if r0.im <= 0.0 && r1.im <= 0.0 { return None; } + if r0.im >= 0.0 && r1.im >= 0.0 { return None; } + + let denom_im = r0.im - r1.im; + if denom_im.abs() < f64::EPSILON { return None; } + let z = r0.im / denom_im; + let x = ((1.0 - z) * r0 + z * r1).re; + if x < ee || x > 1.0 - ee { + return None; + } + + let hit = Complex64::new(p0.0, p0.1) * (1.0 - z) + Complex64::new(q0.0, q0.1) * 0.0; + let _ = hit; + let result = Complex64::new(q0.0, q0.1) * (1.0 - z) + Complex64::new(q1.0, q1.1) * z; + Some((result.re, result.im)) +} + +/// intersection points of an open line segment with an open arc, returning 0, 1, or 2 hits. +pub fn line_arc_intersection( + p0: (f64, f64), + p1: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> Vec<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let a0c = Complex64::new(a0.0, a0.1); + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let c = Complex64::new(cx, cy); + + let d = (p1c - p0c).norm(); + if d < f64::EPSILON { + return Vec::new(); + } + let t = (p1c - p0c) / d; + let v = (c - p0c) / t; + + if v.im.abs() > r { + return Vec::new(); + } + let l = (r * r - v.im * v.im).sqrt(); + let tta = sweep_deg.to_radians().abs(); + + let mut out = Vec::with_capacity(2); + + if (l / r) < 1.0e-5 { + let hit = p0c + Complex64::new(v.re, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + return out; + } + + for sign in [1.0, -1.0] { + let hit = p0c + Complex64::new(v.re + sign * l, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + } + out +} + +/// intersection points of two open arcs, returning 0, 1, or 2 hits. +pub fn arc_arc_intersection( + a0_p0: (f64, f64), + a0_p1: (f64, f64), + a0_sweep_deg: f64, + a1_p0: (f64, f64), + a1_p1: (f64, f64), + a1_sweep_deg: f64, +) -> Vec<(f64, f64)> { + let (c0x, c0y, r0) = circle_from_arc(a0_p0, a0_p1, a0_sweep_deg); + let (c1x, c1y, r1) = circle_from_arc(a1_p0, a1_p1, a1_sweep_deg); + let c0 = Complex64::new(c0x, c0y); + let c1 = Complex64::new(c1x, c1y); + let d = (c1 - c0).norm(); + + if d > r0 + r1 || d < 1.0e-8 { + return Vec::new(); + } + + let l = ((r0 + r1 - d) * (d + r0 - r1) * (d - r0 + r1) * (d + r0 + r1)).sqrt() / (2.0 * d); + let c = 1.0 + (r0 / d) * (r0 / d) - (r1 / d) * (r1 / d); + let t = (c1 - c0) / d; + let tta0 = a0_sweep_deg.to_radians().abs(); + let tta1 = a1_sweep_deg.to_radians().abs(); + let a0c = Complex64::new(a0_p0.0, a0_p0.1); + let a1c = Complex64::new(a1_p0.0, a1_p0.1); + + let mut out = Vec::with_capacity(2); + + let first = c0 + Complex64::new(c * d / 2.0, l) * t; + let z0 = ((first - c0) / (a0c - c0)).arg(); + let z1 = ((first - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((first.re, first.im)); + } + + // tangent-touch case: only one intersection + if (d - r0 + r1).abs() / (r0 + r1) < 1.0e-5 { + return out; + } + + let second = c0 + Complex64::new(c * d / 2.0, -l) * t; + let z0 = ((second - c0) / (a0c - c0)).arg(); + let z1 = ((second - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((second.re, second.im)); + } + + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: (f64, f64), b: (f64, f64), tol: f64) -> bool { + ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt() < tol + } + + #[test] + fn line_line_crosses_at_origin() { + // (-1,-1)->(1,1) crosses (-1,1)->(1,-1) at origin. + let hit = line_line_intersection((-1.0, -1.0), (1.0, 1.0), (-1.0, 1.0), (1.0, -1.0)); + assert!(hit.is_some()); + assert!(close(hit.unwrap(), (0.0, 0.0), 1e-9)); + } + + #[test] + fn line_line_parallel_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn line_line_t_junction_at_endpoint_rejects() { + // the second segment's q0 sits ON the first segment's interior + // but the prospective endpoint check should reject (intersection at q0). + let hit = line_line_intersection((0.0, 0.0), (2.0, 0.0), (1.0, 0.0), (1.0, 1.0)); + // shape: one endpoint sits on the other segment, no clean crossing. + // GetIntersection returns FALSE when an intersection lies within ee of an endpoint of the + // prospective line — here the meeting point IS the prospective q0, so no split is wanted. + assert!(hit.is_none()); + } + + #[test] + fn line_line_shared_endpoint_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (1.0, 0.0), (2.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn circle_from_quarter_arc_unit() { + // a quarter circle (90 deg sweep) from (1, 0) to (0, 1) sits on the unit circle centered at origin. + let (cx, cy, r) = circle_from_arc((1.0, 0.0), (0.0, 1.0), 90.0); + assert!((cx - 0.0).abs() < 1e-9); + assert!((cy - 0.0).abs() < 1e-9); + assert!((r - 1.0).abs() < 1e-9); + } + + #[test] + fn line_through_quarter_arc_crosses_once() { + // diagonal line from (-2, 0.5) to (2, 0.5) crosses the quarter circle (1,0)->(0,1) at one point. + let hits = line_arc_intersection((-2.0, 0.5), (2.0, 0.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert_eq!(hits.len(), 1); + let h = hits[0]; + assert!((h.0 * h.0 + h.1 * h.1 - 1.0).abs() < 1e-9, "point should lie on unit circle: {h:?}"); + assert!((h.1 - 0.5).abs() < 1e-9); + assert!(h.0 > 0.0); + } + + #[test] + fn line_outside_arc_does_not_cross() { + // line clearly outside the arc swept region. + let hits = line_arc_intersection((-2.0, 2.5), (2.0, 2.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!(hits.is_empty()); + } + + #[test] + fn arc_arc_two_intersections() { + // two unit circles offset by 1.0 along x, both swept 360 degrees... but we use partial arcs + // shaped to clip both crossings in their swept range. + // circle A: center (0,0), arcs covering full top half (180 deg from (1,0) to (-1,0)). + // circle B: center (1,0), arcs covering full top half (180 deg from (2,0) to (0,0)). + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (2.0, 0.0), (0.0, 0.0), 180.0, + ); + assert_eq!(hits.len(), 1, "two top-half arcs of overlapping unit circles meet at one upper crossing: {hits:?}"); + let h = hits[0]; + assert!((h.0 - 0.5).abs() < 1e-9, "x = 0.5 by symmetry, got {h:?}"); + assert!(h.1 > 0.0, "intersection should be above the x-axis: {h:?}"); + } + + #[test] + fn arc_arc_disjoint() { + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (11.0, 0.0), (9.0, 0.0), 180.0, + ); + assert!(hits.is_empty()); + } + + #[test] + fn shortest_distance_from_segment_perpendicular() { + // point (0, 1) to segment (-1, 0)->(1, 0) is distance 1. + let d = shortest_distance_from_segment((0.0, 1.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_segment_clamps_to_endpoint() { + // point (-2, 0) to segment (-1, 0)->(1, 0) is distance 1 (clamps to left endpoint). + let d = shortest_distance_from_segment((-2.0, 0.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_arc_radial() { + // point at origin distance to unit-circle quarter arc is 1. + let d = shortest_distance_from_arc((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!((d - 1.0).abs() < 1e-9); + } +} diff --git a/crates/femm-doc-curr/src/lib.rs b/crates/femm-doc-curr/src/lib.rs index cf6f41f..8f61f72 100644 --- a/crates/femm-doc-curr/src/lib.rs +++ b/crates/femm-doc-curr/src/lib.rs @@ -1,9 +1,11 @@ //! current-flow pre-processor document model: geometry, properties, parser, writer. pub mod geom; +pub mod geom_math; pub mod props; pub mod parser; pub mod writer; +pub mod edit; use num_complex::Complex64; diff --git a/crates/femm-doc-elec/src/edit.rs b/crates/femm-doc-elec/src/edit.rs new file mode 100644 index 0000000..7e87bbd --- /dev/null +++ b/crates/femm-doc-elec/src/edit.rs @@ -0,0 +1,612 @@ +//! geometry editing primitives on [`FemmDoc`]: add, delete, closest-point queries. + +use crate::geom_math::{ + arc_arc_intersection, circle_from_arc, line_arc_intersection, line_line_intersection, + shortest_distance_from_arc, shortest_distance_from_segment, +}; +use crate::{ArcSegment, BlockLabel, FemmDoc, Node, Segment}; +use num_complex::Complex64; + +/// fraction of the node-bbox diagonal used as auto-tolerance for intersection-node coalescing. +const BBOX_TOLERANCE_FRAC: f64 = 1.0e-6; +/// fraction of a segment's length within which an off-endpoint node triggers a recursive split. +const ON_LINE_FRAC: f64 = 1.0e-5; + +impl FemmDoc { + /// adds a node at (x, y), returning the index of an existing node within `tol` distance when present. + pub fn add_node(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, n) in self.nodes.iter().enumerate() { + let dx = n.x - x; + let dy = n.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.nodes.push(Node { + x, y, + boundary_marker: String::new(), + in_conductor: String::new(), + in_group: 0, + selected: false, + }); + self.nodes.len() - 1 + } + + /// adds a block label at (x, y), returning the index of an existing label within `tol` when present. + pub fn add_block_label(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, b) in self.block_labels.iter().enumerate() { + let dx = b.x - x; + let dy = b.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.block_labels.push(BlockLabel { + x, y, + max_area: 0.0, + block_type: String::from(""), + in_conductor: String::from(""), + in_group: 0, + is_external: false, + is_default: false, + selected: false, + }); + self.block_labels.len() - 1 + } + + /// adds a segment between two node indices, splitting at every crossing and through any on-line node. + pub fn add_segment(&mut self, n0: i32, n1: i32) -> bool { + self.add_segment_with_marker(n0, n1, "") + } + + /// PSLG-aware variant propagating a boundary-marker name onto every resulting segment piece. + pub fn add_segment_with_marker(&mut self, n0: i32, n1: i32, marker: &str) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // duplicate either-orientation. + let (a, b) = if n0 < n1 { (n0, n1) } else { (n1, n0) }; + for s in &self.segments { + let (sa, sb) = if s.n0 < s.n1 { (s.n0, s.n1) } else { (s.n1, s.n0) }; + if sa == a && sb == b { return false; } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + if let Some(hit) = line_line_intersection(n0p, n1p, sp0, sp1) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in line_arc_intersection(n0p, n1p, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + self.segments.push(Segment { + n0, n1, + max_side_length: -1.0, + boundary_marker: marker.to_string(), + in_conductor: String::new(), + hidden: false, + in_group: 0, + selected: false, + }); + + // first non-endpoint node on the new segment's interior, if any. + let length = ((n1p.0 - n0p.0).powi(2) + (n1p.1 - n0p.1).powi(2)).sqrt(); + let dmin = length * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_segment(np, n0p, n1p); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.segments.pop(); + let a = self.add_segment_with_marker(n0, mid, marker); + let b = self.add_segment_with_marker(mid, n1, marker); + a || b + } else { + true + } + } + + /// auto-tolerance for intersection-node coalescing, derived from the node bounding box. + fn bbox_tolerance(&self) -> f64 { + nodes_bbox_tolerance(&self.nodes) + } + + /// rebuilds every list through the PSLG-aware add primitives, catching crossings missed by incremental edits. + pub fn enforce_pslg(&mut self) { + let old_nodes = std::mem::take(&mut self.nodes); + let old_segments = std::mem::take(&mut self.segments); + let old_arcs = std::mem::take(&mut self.arcs); + let old_block_labels = std::mem::take(&mut self.block_labels); + + let tol = nodes_bbox_tolerance(&old_nodes); + + // dedupes by position with metadata preserved. + for n in &old_nodes { + let mut duplicate = false; + for existing in &self.nodes { + if (existing.x - n.x).hypot(existing.y - n.y) < tol { + duplicate = true; + break; + } + } + if !duplicate { + self.nodes.push(n.clone()); + } + } + + for s in old_segments { + let p0 = (old_nodes[s.n0 as usize].x, old_nodes[s.n0 as usize].y); + let p1 = (old_nodes[s.n1 as usize].x, old_nodes[s.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_segment_with_marker(n0 as i32, n1 as i32, &s.boundary_marker); + } + + for a in old_arcs { + let p0 = (old_nodes[a.n0 as usize].x, old_nodes[a.n0 as usize].y); + let p1 = (old_nodes[a.n1 as usize].x, old_nodes[a.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_arc_segment_with_template(n0 as i32, n1 as i32, a.arc_length, &a); + } + + self.block_labels = old_block_labels; + } + + /// adds an arc between two node indices, splitting at every crossing and through any on-arc node. + pub fn add_arc_segment(&mut self, n0: i32, n1: i32, arc_length_deg: f64) -> bool { + let template = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + max_side_length: 10.0, + boundary_marker: String::new(), + in_conductor: String::new(), + hidden: false, + in_group: 0, + normal_direction: true, + selected: false, + }; + self.add_arc_segment_with_template(n0, n1, arc_length_deg, &template) + } + + /// PSLG-aware variant propagating boundary marker and side-length metadata onto every arc piece. + pub fn add_arc_segment_with_template( + &mut self, + n0: i32, + n1: i32, + arc_length_deg: f64, + template: &ArcSegment, + ) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // same directed endpoints with similar sweep counts as duplicate. + for a in &self.arcs { + if a.n0 == n0 && a.n1 == n1 && (a.arc_length - arc_length_deg).abs() < 1.0e-2 { + return false; + } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + for hit in line_arc_intersection(sp0, sp1, n0p, n1p, arc_length_deg) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in arc_arc_intersection(n0p, n1p, arc_length_deg, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + let new_arc = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + ..template.clone() + }; + self.arcs.push(new_arc); + + // first non-endpoint node on the new arc's sweep range, if any. + let (cx, cy, radius) = circle_from_arc(n0p, n1p, arc_length_deg); + let sweep_rad = arc_length_deg.to_radians(); + let arc_length_world = radius * sweep_rad.abs(); + let dmin = arc_length_world * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_arc(np, n0p, n1p, arc_length_deg); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.arcs.pop(); + let mid_pos = (self.nodes[mid as usize].x, self.nodes[mid as usize].y); + let c = Complex64::new(cx, cy); + let a0 = Complex64::new(n0p.0, n0p.1); + let a1 = Complex64::new(n1p.0, n1p.1); + let a2 = Complex64::new(mid_pos.0, mid_pos.1); + let sweep_to_mid = ((a2 - c) / (a0 - c)).arg().to_degrees(); + let sweep_from_mid = ((a1 - c) / (a2 - c)).arg().to_degrees(); + let a = self.add_arc_segment_with_template(n0, mid, sweep_to_mid, template); + let b = self.add_arc_segment_with_template(mid, n1, sweep_from_mid, template); + a || b + } else { + true + } + } + + /// removes selected nodes and rewrites segment/arc endpoint indices to drop references. + pub fn delete_selected_nodes(&mut self) -> usize { + let keep: Vec = self.nodes.iter().map(|n| !n.selected).collect(); + let mut remap: Vec = Vec::with_capacity(keep.len()); + let mut next: i32 = 0; + for &k in &keep { + remap.push(if k { let r = next; next += 1; r } else { -1 }); + } + let removed = self.nodes.len() - next as usize; + + let mut new_nodes = Vec::with_capacity(next as usize); + for (i, n) in self.nodes.drain(..).enumerate() { + if keep[i] { new_nodes.push(n); } + } + self.nodes = new_nodes; + + self.segments.retain(|s| { + let a = s.n0 as usize; + let b = s.n1 as usize; + a < keep.len() && b < keep.len() && keep[a] && keep[b] + }); + for s in &mut self.segments { + s.n0 = remap[s.n0 as usize]; + s.n1 = remap[s.n1 as usize]; + } + + self.arcs.retain(|a| { + let i = a.n0 as usize; + let j = a.n1 as usize; + i < keep.len() && j < keep.len() && keep[i] && keep[j] + }); + for a in &mut self.arcs { + a.n0 = remap[a.n0 as usize]; + a.n1 = remap[a.n1 as usize]; + } + + removed + } + + /// removes selected segments and returns the count. + pub fn delete_selected_segments(&mut self) -> usize { + let before = self.segments.len(); + self.segments.retain(|s| !s.selected); + before - self.segments.len() + } + + /// removes selected arc segments and returns the count. + pub fn delete_selected_arcs(&mut self) -> usize { + let before = self.arcs.len(); + self.arcs.retain(|a| !a.selected); + before - self.arcs.len() + } + + /// removes selected block labels and returns the count. + pub fn delete_selected_block_labels(&mut self) -> usize { + let before = self.block_labels.len(); + self.block_labels.retain(|b| !b.selected); + before - self.block_labels.len() + } + + /// returns the index of the closest node to (x, y), or None when the node list is empty. + pub fn closest_node(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, n) in self.nodes.iter().enumerate() { + let d = (n.x - x).hypot(n.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the closest block label to (x, y), or None when none exist. + pub fn closest_block_label(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, b) in self.block_labels.iter().enumerate() { + let d = (b.x - x).hypot(b.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the segment whose nearest point is closest to (x, y). + pub fn closest_segment(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, s) in self.segments.iter().enumerate() { + let (Some(p0), Some(p1)) = (self.nodes.get(s.n0 as usize), self.nodes.get(s.n1 as usize)) else { continue }; + let d = point_to_segment_distance(x, y, p0.x, p0.y, p1.x, p1.y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// clears the selection flag on every geometric entity in the doc. + pub fn clear_selection(&mut self) { + for n in &mut self.nodes { n.selected = false; } + for s in &mut self.segments { s.selected = false; } + for a in &mut self.arcs { a.selected = false; } + for b in &mut self.block_labels { b.selected = false; } + } +} + +/// auto-tolerance derived from the bounding box of a node slice. +fn nodes_bbox_tolerance(nodes: &[Node]) -> f64 { + if nodes.len() < 2 { + return 1.0e-8; + } + let mut xmin = f64::INFINITY; + let mut xmax = f64::NEG_INFINITY; + let mut ymin = f64::INFINITY; + let mut ymax = f64::NEG_INFINITY; + for n in nodes { + xmin = xmin.min(n.x); xmax = xmax.max(n.x); + ymin = ymin.min(n.y); ymax = ymax.max(n.y); + } + let dx = xmax - xmin; + let dy = ymax - ymin; + (dx * dx + dy * dy).sqrt() * BBOX_TOLERANCE_FRAC +} + +/// Euclidean distance from (px, py) to the segment between (ax, ay) and (bx, by). +fn point_to_segment_distance(px: f64, py: f64, ax: f64, ay: f64, bx: f64, by: f64) -> f64 { + let dx = bx - ax; + let dy = by - ay; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return (px - ax).hypot(py - ay); + } + let t = (((px - ax) * dx + (py - ay) * dy) / len2).clamp(0.0, 1.0); + let cx = ax + t * dx; + let cy = ay + t * dy; + (px - cx).hypot(py - cy) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn doc_with_corners() -> FemmDoc { + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0,-1.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d + } + + #[test] + fn duplicate_segment_rejected() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert!(!d.add_segment(0, 1)); + assert!(!d.add_segment(1, 0)); + assert_eq!(d.segments.len(), 1); + } + + #[test] + fn segment_passing_through_existing_node_splits() { + // three colinear nodes (-1,0), (1,0), (0,0). adding the outer-to-outer segment + // splits at the midpoint node, producing two pieces meeting there. + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert_eq!(d.segments.len(), 2); + let touches_mid = d.segments.iter().filter(|s| s.n0 == 2 || s.n1 == 2).count(); + assert_eq!(touches_mid, 2); + } + + #[test] + fn enforce_pslg_splits_first_segment_at_late_intersection() { + // incremental edit splits only the second of two crossing segments. enforce_pslg + // back-splits the first segment through the intersection node. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + assert_eq!(d.segments.len(), 3, "pre-enforce: one whole, one split"); + + d.enforce_pslg(); + + assert_eq!(d.nodes.len(), 5, "no node added or merged by enforce"); + assert_eq!(d.segments.len(), 4, "both segments split at the origin"); + let touches_origin = d.segments.iter().filter(|s| s.n0 == 4 || s.n1 == 4).count(); + assert_eq!(touches_origin, 4); + } + + #[test] + fn enforce_pslg_preserves_segment_boundary_marker() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment_with_marker(0, 1, "outer")); + + d.enforce_pslg(); + + assert_eq!(d.segments.len(), 1); + assert_eq!(d.segments[0].boundary_marker, "outer"); + } + + #[test] + fn enforce_pslg_preserves_arc_metadata() { + // arc carrying non-default marker, side length, group, and direction. + // enforce_pslg round-trips every field intact. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let template = ArcSegment { + n0: 0, + n1: 1, + arc_length: 90.0, + max_side_length: 5.0, + boundary_marker: String::from("outer"), + in_conductor: String::new(), + hidden: true, + in_group: 7, + normal_direction: false, + selected: false, + }; + assert!(d.add_arc_segment_with_template(0, 1, 90.0, &template)); + + d.enforce_pslg(); + + assert_eq!(d.arcs.len(), 1); + let a = &d.arcs[0]; + assert_eq!(a.boundary_marker, "outer"); + assert!((a.max_side_length - 5.0).abs() < 1e-12); + assert!(a.hidden); + assert_eq!(a.in_group, 7); + assert!(!a.normal_direction); + } + + #[test] + fn duplicate_arc_rejected() { + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(!d.add_arc_segment(0, 1, 90.0)); + // reversed endpoints define a distinct arc. + assert!(d.add_arc_segment(1, 0, 90.0)); + assert_eq!(d.arcs.len(), 2); + } + + #[test] + fn add_segment_crossing_unit_quarter_arc_splits_segment() { + // unit quarter arc from (1,0) to (0,1); horizontal line at y=0.5 crosses the arc once. + // expect a fresh intersection node near (0.866, 0.5) and the segment split into two. + let mut d = FemmDoc::default(); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d.add_node(-2.0, 0.5, 0.0); + d.add_node( 2.0, 0.5, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5); + let n4 = &d.nodes[4]; + assert!((n4.x - 0.75_f64.sqrt()).abs() < 1e-6); + assert!((n4.y - 0.5).abs() < 1e-9); + assert_eq!(d.segments.len(), 2); + assert_eq!(d.arcs.len(), 1); + } + + #[test] + fn add_arc_passing_through_existing_node_splits_into_two_arcs() { + // pre-existing node at (cos45, sin45) lies on the unit quarter arc from (1,0) to (0,1). + // adding the arc detects the on-arc node and emits two sub-arcs of 45 degrees each. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let mid_x = std::f64::consts::FRAC_1_SQRT_2; + d.add_node(mid_x, mid_x, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + + assert_eq!(d.arcs.len(), 2); + let touches_mid = d.arcs.iter().filter(|a| a.n0 == 2 || a.n1 == 2).count(); + assert_eq!(touches_mid, 2); + // split sweeps sum to the parent 90 degrees, within floating slop. + let total: f64 = d.arcs.iter().map(|a| a.arc_length).sum(); + assert!((total - 90.0).abs() < 1e-6, "sum of split sweeps = {total}"); + } + + #[test] + fn crossing_second_segment_splits_at_intersection() { + // two perpendicular segments through the origin. horizontal first, vertical second. + // the vertical add inserts an origin node and self-splits; the horizontal stays whole. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5, "intersection node added"); + let new_idx = 4; + assert!((d.nodes[new_idx].x).abs() < 1e-9); + assert!((d.nodes[new_idx].y).abs() < 1e-9); + + // expected pieces: 0->1 (horizontal whole), 2->4 and 4->3 (vertical halves). + assert_eq!(d.segments.len(), 3); + let touches_new = d.segments.iter().filter(|s| s.n0 == new_idx as i32 || s.n1 == new_idx as i32).count(); + assert_eq!(touches_new, 2); + } +} diff --git a/crates/femm-doc-elec/src/geom_math.rs b/crates/femm-doc-elec/src/geom_math.rs new file mode 100644 index 0000000..37eeb54 --- /dev/null +++ b/crates/femm-doc-elec/src/geom_math.rs @@ -0,0 +1,322 @@ +//! 2D intersection math and distance queries over (x, y) pairs and degree-sweep arcs. + +use num_complex::Complex64; + +/// closeness factor for common-endpoint and small-interval rejection. +const EPS_FRAC: f64 = 1.0e-8; + +/// recovers the center and radius of the circle through an arc's endpoints with the given sweep. +pub fn circle_from_arc(p0: (f64, f64), p1: (f64, f64), sweep_deg: f64) -> (f64, f64, f64) { + let a0 = Complex64::new(p0.0, p0.1); + let a1 = Complex64::new(p1.0, p1.1); + let chord = (a1 - a0).norm(); + let t = (a1 - a0) / chord; + let tta = sweep_deg.to_radians(); + let r = chord / (2.0 * (tta / 2.0).sin()); + let c = a0 + (Complex64::new(chord / 2.0, (r * r - chord * chord / 4.0).sqrt())) * t; + (c.re, c.im, r.abs()) +} + +/// shortest distance from a point to a segment. +pub fn shortest_distance_from_segment(p: (f64, f64), a: (f64, f64), b: (f64, f64)) -> f64 { + let dx = b.0 - a.0; + let dy = b.1 - a.1; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return ((p.0 - a.0).powi(2) + (p.1 - a.1).powi(2)).sqrt(); + } + let t = (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / len2).clamp(0.0, 1.0); + let cx = a.0 + t * dx; + let cy = a.1 + t * dy; + ((p.0 - cx).powi(2) + (p.1 - cy).powi(2)).sqrt() +} + +/// shortest distance from a point to a sweep-defined arc. +pub fn shortest_distance_from_arc( + p: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> f64 { + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let pc = Complex64::new(p.0 - cx, p.1 - cy); + let d = pc.norm(); + if d == 0.0 { + return r; + } + let t = pc / d; + let foot = pc - Complex64::new(r * t.re, r * t.im); + let l = foot.norm(); + let a0_dir = Complex64::new(a0.0 - cx, a0.1 - cy); + let z = ((t / a0_dir).arg() * 180.0 / std::f64::consts::PI + 360.0) % 360.0; + let sweep_abs = sweep_deg.abs(); + if z > 0.0 && z < sweep_abs { + return l; + } + let e0 = ((p.0 - a0.0).powi(2) + (p.1 - a0.1).powi(2)).sqrt(); + let e1 = ((p.0 - a1.0).powi(2) + (p.1 - a1.1).powi(2)).sqrt(); + e0.min(e1) +} + +/// intersection of the open segments p0->p1 and q0->q1, or None when none exists. +pub fn line_line_intersection( + p0: (f64, f64), + p1: (f64, f64), + q0: (f64, f64), + q1: (f64, f64), +) -> Option<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let q0c = Complex64::new(q0.0, q0.1); + let q1c = Complex64::new(q1.0, q1.1); + + // shared endpoint -> no other intersection + if (p0c - q0c).norm() < f64::EPSILON + || (p0c - q1c).norm() < f64::EPSILON + || (p1c - q0c).norm() < f64::EPSILON + || (p1c - q1c).norm() < f64::EPSILON + { + return None; + } + + let ee = ((p1c - p0c).norm()).min((q1c - q0c).norm()) * EPS_FRAC; + + let denom = p1c - p0c; + if denom.norm() < f64::EPSILON { + return None; + } + let r0 = (q0c - p0c) / denom; + let r1 = (q1c - p0c) / denom; + + if r0.re <= 0.0 && r1.re <= 0.0 { return None; } + if r0.re >= 1.0 && r1.re >= 1.0 { return None; } + if r0.im <= 0.0 && r1.im <= 0.0 { return None; } + if r0.im >= 0.0 && r1.im >= 0.0 { return None; } + + let denom_im = r0.im - r1.im; + if denom_im.abs() < f64::EPSILON { return None; } + let z = r0.im / denom_im; + let x = ((1.0 - z) * r0 + z * r1).re; + if x < ee || x > 1.0 - ee { + return None; + } + + let hit = Complex64::new(p0.0, p0.1) * (1.0 - z) + Complex64::new(q0.0, q0.1) * 0.0; + let _ = hit; + let result = Complex64::new(q0.0, q0.1) * (1.0 - z) + Complex64::new(q1.0, q1.1) * z; + Some((result.re, result.im)) +} + +/// intersection points of an open line segment with an open arc, returning 0, 1, or 2 hits. +pub fn line_arc_intersection( + p0: (f64, f64), + p1: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> Vec<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let a0c = Complex64::new(a0.0, a0.1); + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let c = Complex64::new(cx, cy); + + let d = (p1c - p0c).norm(); + if d < f64::EPSILON { + return Vec::new(); + } + let t = (p1c - p0c) / d; + let v = (c - p0c) / t; + + if v.im.abs() > r { + return Vec::new(); + } + let l = (r * r - v.im * v.im).sqrt(); + let tta = sweep_deg.to_radians().abs(); + + let mut out = Vec::with_capacity(2); + + if (l / r) < 1.0e-5 { + let hit = p0c + Complex64::new(v.re, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + return out; + } + + for sign in [1.0, -1.0] { + let hit = p0c + Complex64::new(v.re + sign * l, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + } + out +} + +/// intersection points of two open arcs, returning 0, 1, or 2 hits. +pub fn arc_arc_intersection( + a0_p0: (f64, f64), + a0_p1: (f64, f64), + a0_sweep_deg: f64, + a1_p0: (f64, f64), + a1_p1: (f64, f64), + a1_sweep_deg: f64, +) -> Vec<(f64, f64)> { + let (c0x, c0y, r0) = circle_from_arc(a0_p0, a0_p1, a0_sweep_deg); + let (c1x, c1y, r1) = circle_from_arc(a1_p0, a1_p1, a1_sweep_deg); + let c0 = Complex64::new(c0x, c0y); + let c1 = Complex64::new(c1x, c1y); + let d = (c1 - c0).norm(); + + if d > r0 + r1 || d < 1.0e-8 { + return Vec::new(); + } + + let l = ((r0 + r1 - d) * (d + r0 - r1) * (d - r0 + r1) * (d + r0 + r1)).sqrt() / (2.0 * d); + let c = 1.0 + (r0 / d) * (r0 / d) - (r1 / d) * (r1 / d); + let t = (c1 - c0) / d; + let tta0 = a0_sweep_deg.to_radians().abs(); + let tta1 = a1_sweep_deg.to_radians().abs(); + let a0c = Complex64::new(a0_p0.0, a0_p0.1); + let a1c = Complex64::new(a1_p0.0, a1_p0.1); + + let mut out = Vec::with_capacity(2); + + let first = c0 + Complex64::new(c * d / 2.0, l) * t; + let z0 = ((first - c0) / (a0c - c0)).arg(); + let z1 = ((first - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((first.re, first.im)); + } + + // tangent-touch case: only one intersection + if (d - r0 + r1).abs() / (r0 + r1) < 1.0e-5 { + return out; + } + + let second = c0 + Complex64::new(c * d / 2.0, -l) * t; + let z0 = ((second - c0) / (a0c - c0)).arg(); + let z1 = ((second - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((second.re, second.im)); + } + + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: (f64, f64), b: (f64, f64), tol: f64) -> bool { + ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt() < tol + } + + #[test] + fn line_line_crosses_at_origin() { + // (-1,-1)->(1,1) crosses (-1,1)->(1,-1) at origin. + let hit = line_line_intersection((-1.0, -1.0), (1.0, 1.0), (-1.0, 1.0), (1.0, -1.0)); + assert!(hit.is_some()); + assert!(close(hit.unwrap(), (0.0, 0.0), 1e-9)); + } + + #[test] + fn line_line_parallel_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn line_line_t_junction_at_endpoint_rejects() { + // the second segment's q0 sits ON the first segment's interior + // but the prospective endpoint check should reject (intersection at q0). + let hit = line_line_intersection((0.0, 0.0), (2.0, 0.0), (1.0, 0.0), (1.0, 1.0)); + // shape: one endpoint sits on the other segment, no clean crossing. + // GetIntersection returns FALSE when an intersection lies within ee of an endpoint of the + // prospective line — here the meeting point IS the prospective q0, so no split is wanted. + assert!(hit.is_none()); + } + + #[test] + fn line_line_shared_endpoint_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (1.0, 0.0), (2.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn circle_from_quarter_arc_unit() { + // a quarter circle (90 deg sweep) from (1, 0) to (0, 1) sits on the unit circle centered at origin. + let (cx, cy, r) = circle_from_arc((1.0, 0.0), (0.0, 1.0), 90.0); + assert!((cx - 0.0).abs() < 1e-9); + assert!((cy - 0.0).abs() < 1e-9); + assert!((r - 1.0).abs() < 1e-9); + } + + #[test] + fn line_through_quarter_arc_crosses_once() { + // diagonal line from (-2, 0.5) to (2, 0.5) crosses the quarter circle (1,0)->(0,1) at one point. + let hits = line_arc_intersection((-2.0, 0.5), (2.0, 0.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert_eq!(hits.len(), 1); + let h = hits[0]; + assert!((h.0 * h.0 + h.1 * h.1 - 1.0).abs() < 1e-9, "point should lie on unit circle: {h:?}"); + assert!((h.1 - 0.5).abs() < 1e-9); + assert!(h.0 > 0.0); + } + + #[test] + fn line_outside_arc_does_not_cross() { + // line clearly outside the arc swept region. + let hits = line_arc_intersection((-2.0, 2.5), (2.0, 2.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!(hits.is_empty()); + } + + #[test] + fn arc_arc_two_intersections() { + // two unit circles offset by 1.0 along x, both swept 360 degrees... but we use partial arcs + // shaped to clip both crossings in their swept range. + // circle A: center (0,0), arcs covering full top half (180 deg from (1,0) to (-1,0)). + // circle B: center (1,0), arcs covering full top half (180 deg from (2,0) to (0,0)). + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (2.0, 0.0), (0.0, 0.0), 180.0, + ); + assert_eq!(hits.len(), 1, "two top-half arcs of overlapping unit circles meet at one upper crossing: {hits:?}"); + let h = hits[0]; + assert!((h.0 - 0.5).abs() < 1e-9, "x = 0.5 by symmetry, got {h:?}"); + assert!(h.1 > 0.0, "intersection should be above the x-axis: {h:?}"); + } + + #[test] + fn arc_arc_disjoint() { + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (11.0, 0.0), (9.0, 0.0), 180.0, + ); + assert!(hits.is_empty()); + } + + #[test] + fn shortest_distance_from_segment_perpendicular() { + // point (0, 1) to segment (-1, 0)->(1, 0) is distance 1. + let d = shortest_distance_from_segment((0.0, 1.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_segment_clamps_to_endpoint() { + // point (-2, 0) to segment (-1, 0)->(1, 0) is distance 1 (clamps to left endpoint). + let d = shortest_distance_from_segment((-2.0, 0.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_arc_radial() { + // point at origin distance to unit-circle quarter arc is 1. + let d = shortest_distance_from_arc((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!((d - 1.0).abs() < 1e-9); + } +} diff --git a/crates/femm-doc-elec/src/lib.rs b/crates/femm-doc-elec/src/lib.rs index 95276e4..4e1b98c 100644 --- a/crates/femm-doc-elec/src/lib.rs +++ b/crates/femm-doc-elec/src/lib.rs @@ -1,9 +1,11 @@ //! electrostatic pre-processor document model: geometry, properties, parser, writer. pub mod geom; +pub mod geom_math; pub mod props; pub mod parser; pub mod writer; +pub mod edit; pub use geom::{ArcSegment, BlockLabel, Node, Segment}; pub use props::{BoundaryProp, ConductorProp, MaterialProp, PointProp}; diff --git a/crates/femm-doc-heat/src/edit.rs b/crates/femm-doc-heat/src/edit.rs new file mode 100644 index 0000000..7e87bbd --- /dev/null +++ b/crates/femm-doc-heat/src/edit.rs @@ -0,0 +1,612 @@ +//! geometry editing primitives on [`FemmDoc`]: add, delete, closest-point queries. + +use crate::geom_math::{ + arc_arc_intersection, circle_from_arc, line_arc_intersection, line_line_intersection, + shortest_distance_from_arc, shortest_distance_from_segment, +}; +use crate::{ArcSegment, BlockLabel, FemmDoc, Node, Segment}; +use num_complex::Complex64; + +/// fraction of the node-bbox diagonal used as auto-tolerance for intersection-node coalescing. +const BBOX_TOLERANCE_FRAC: f64 = 1.0e-6; +/// fraction of a segment's length within which an off-endpoint node triggers a recursive split. +const ON_LINE_FRAC: f64 = 1.0e-5; + +impl FemmDoc { + /// adds a node at (x, y), returning the index of an existing node within `tol` distance when present. + pub fn add_node(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, n) in self.nodes.iter().enumerate() { + let dx = n.x - x; + let dy = n.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.nodes.push(Node { + x, y, + boundary_marker: String::new(), + in_conductor: String::new(), + in_group: 0, + selected: false, + }); + self.nodes.len() - 1 + } + + /// adds a block label at (x, y), returning the index of an existing label within `tol` when present. + pub fn add_block_label(&mut self, x: f64, y: f64, tol: f64) -> usize { + if tol > 0.0 { + for (i, b) in self.block_labels.iter().enumerate() { + let dx = b.x - x; + let dy = b.y - y; + if (dx * dx + dy * dy).sqrt() <= tol { + return i; + } + } + } + self.block_labels.push(BlockLabel { + x, y, + max_area: 0.0, + block_type: String::from(""), + in_conductor: String::from(""), + in_group: 0, + is_external: false, + is_default: false, + selected: false, + }); + self.block_labels.len() - 1 + } + + /// adds a segment between two node indices, splitting at every crossing and through any on-line node. + pub fn add_segment(&mut self, n0: i32, n1: i32) -> bool { + self.add_segment_with_marker(n0, n1, "") + } + + /// PSLG-aware variant propagating a boundary-marker name onto every resulting segment piece. + pub fn add_segment_with_marker(&mut self, n0: i32, n1: i32, marker: &str) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // duplicate either-orientation. + let (a, b) = if n0 < n1 { (n0, n1) } else { (n1, n0) }; + for s in &self.segments { + let (sa, sb) = if s.n0 < s.n1 { (s.n0, s.n1) } else { (s.n1, s.n0) }; + if sa == a && sb == b { return false; } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + if let Some(hit) = line_line_intersection(n0p, n1p, sp0, sp1) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in line_arc_intersection(n0p, n1p, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + self.segments.push(Segment { + n0, n1, + max_side_length: -1.0, + boundary_marker: marker.to_string(), + in_conductor: String::new(), + hidden: false, + in_group: 0, + selected: false, + }); + + // first non-endpoint node on the new segment's interior, if any. + let length = ((n1p.0 - n0p.0).powi(2) + (n1p.1 - n0p.1).powi(2)).sqrt(); + let dmin = length * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_segment(np, n0p, n1p); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.segments.pop(); + let a = self.add_segment_with_marker(n0, mid, marker); + let b = self.add_segment_with_marker(mid, n1, marker); + a || b + } else { + true + } + } + + /// auto-tolerance for intersection-node coalescing, derived from the node bounding box. + fn bbox_tolerance(&self) -> f64 { + nodes_bbox_tolerance(&self.nodes) + } + + /// rebuilds every list through the PSLG-aware add primitives, catching crossings missed by incremental edits. + pub fn enforce_pslg(&mut self) { + let old_nodes = std::mem::take(&mut self.nodes); + let old_segments = std::mem::take(&mut self.segments); + let old_arcs = std::mem::take(&mut self.arcs); + let old_block_labels = std::mem::take(&mut self.block_labels); + + let tol = nodes_bbox_tolerance(&old_nodes); + + // dedupes by position with metadata preserved. + for n in &old_nodes { + let mut duplicate = false; + for existing in &self.nodes { + if (existing.x - n.x).hypot(existing.y - n.y) < tol { + duplicate = true; + break; + } + } + if !duplicate { + self.nodes.push(n.clone()); + } + } + + for s in old_segments { + let p0 = (old_nodes[s.n0 as usize].x, old_nodes[s.n0 as usize].y); + let p1 = (old_nodes[s.n1 as usize].x, old_nodes[s.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_segment_with_marker(n0 as i32, n1 as i32, &s.boundary_marker); + } + + for a in old_arcs { + let p0 = (old_nodes[a.n0 as usize].x, old_nodes[a.n0 as usize].y); + let p1 = (old_nodes[a.n1 as usize].x, old_nodes[a.n1 as usize].y); + let (Some(n0), Some(n1)) = ( + self.closest_node(p0.0, p0.1), + self.closest_node(p1.0, p1.1), + ) else { continue }; + self.add_arc_segment_with_template(n0 as i32, n1 as i32, a.arc_length, &a); + } + + self.block_labels = old_block_labels; + } + + /// adds an arc between two node indices, splitting at every crossing and through any on-arc node. + pub fn add_arc_segment(&mut self, n0: i32, n1: i32, arc_length_deg: f64) -> bool { + let template = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + max_side_length: 10.0, + boundary_marker: String::new(), + in_conductor: String::new(), + hidden: false, + in_group: 0, + normal_direction: true, + selected: false, + }; + self.add_arc_segment_with_template(n0, n1, arc_length_deg, &template) + } + + /// PSLG-aware variant propagating boundary marker and side-length metadata onto every arc piece. + pub fn add_arc_segment_with_template( + &mut self, + n0: i32, + n1: i32, + arc_length_deg: f64, + template: &ArcSegment, + ) -> bool { + if n0 == n1 { return false; } + let nn = self.nodes.len() as i32; + if n0 < 0 || n1 < 0 || n0 >= nn || n1 >= nn { return false; } + + // same directed endpoints with similar sweep counts as duplicate. + for a in &self.arcs { + if a.n0 == n0 && a.n1 == n1 && (a.arc_length - arc_length_deg).abs() < 1.0e-2 { + return false; + } + } + + let n0p = (self.nodes[n0 as usize].x, self.nodes[n0 as usize].y); + let n1p = (self.nodes[n1 as usize].x, self.nodes[n1 as usize].y); + + // crossings with existing segments and arcs. + let mut new_points: Vec<(f64, f64)> = Vec::new(); + for s in &self.segments { + let sp0 = (self.nodes[s.n0 as usize].x, self.nodes[s.n0 as usize].y); + let sp1 = (self.nodes[s.n1 as usize].x, self.nodes[s.n1 as usize].y); + for hit in line_arc_intersection(sp0, sp1, n0p, n1p, arc_length_deg) { + new_points.push(hit); + } + } + for arc in &self.arcs { + let ap0 = (self.nodes[arc.n0 as usize].x, self.nodes[arc.n0 as usize].y); + let ap1 = (self.nodes[arc.n1 as usize].x, self.nodes[arc.n1 as usize].y); + for hit in arc_arc_intersection(n0p, n1p, arc_length_deg, ap0, ap1, arc.arc_length) { + new_points.push(hit); + } + } + + let tol = self.bbox_tolerance(); + for (x, y) in new_points { + self.add_node(x, y, tol); + } + + let new_arc = ArcSegment { + n0, n1, + arc_length: arc_length_deg, + ..template.clone() + }; + self.arcs.push(new_arc); + + // first non-endpoint node on the new arc's sweep range, if any. + let (cx, cy, radius) = circle_from_arc(n0p, n1p, arc_length_deg); + let sweep_rad = arc_length_deg.to_radians(); + let arc_length_world = radius * sweep_rad.abs(); + let dmin = arc_length_world * ON_LINE_FRAC; + + let mut split_at: Option = None; + for (i, node) in self.nodes.iter().enumerate() { + let idx = i as i32; + if idx == n0 || idx == n1 { continue; } + let np = (node.x, node.y); + let de0 = ((np.0 - n0p.0).powi(2) + (np.1 - n0p.1).powi(2)).sqrt(); + let de1 = ((np.0 - n1p.0).powi(2) + (np.1 - n1p.1).powi(2)).sqrt(); + if de0 < dmin || de1 < dmin { continue; } + let d = shortest_distance_from_arc(np, n0p, n1p, arc_length_deg); + if d < dmin { + split_at = Some(idx); + break; + } + } + + if let Some(mid) = split_at { + self.arcs.pop(); + let mid_pos = (self.nodes[mid as usize].x, self.nodes[mid as usize].y); + let c = Complex64::new(cx, cy); + let a0 = Complex64::new(n0p.0, n0p.1); + let a1 = Complex64::new(n1p.0, n1p.1); + let a2 = Complex64::new(mid_pos.0, mid_pos.1); + let sweep_to_mid = ((a2 - c) / (a0 - c)).arg().to_degrees(); + let sweep_from_mid = ((a1 - c) / (a2 - c)).arg().to_degrees(); + let a = self.add_arc_segment_with_template(n0, mid, sweep_to_mid, template); + let b = self.add_arc_segment_with_template(mid, n1, sweep_from_mid, template); + a || b + } else { + true + } + } + + /// removes selected nodes and rewrites segment/arc endpoint indices to drop references. + pub fn delete_selected_nodes(&mut self) -> usize { + let keep: Vec = self.nodes.iter().map(|n| !n.selected).collect(); + let mut remap: Vec = Vec::with_capacity(keep.len()); + let mut next: i32 = 0; + for &k in &keep { + remap.push(if k { let r = next; next += 1; r } else { -1 }); + } + let removed = self.nodes.len() - next as usize; + + let mut new_nodes = Vec::with_capacity(next as usize); + for (i, n) in self.nodes.drain(..).enumerate() { + if keep[i] { new_nodes.push(n); } + } + self.nodes = new_nodes; + + self.segments.retain(|s| { + let a = s.n0 as usize; + let b = s.n1 as usize; + a < keep.len() && b < keep.len() && keep[a] && keep[b] + }); + for s in &mut self.segments { + s.n0 = remap[s.n0 as usize]; + s.n1 = remap[s.n1 as usize]; + } + + self.arcs.retain(|a| { + let i = a.n0 as usize; + let j = a.n1 as usize; + i < keep.len() && j < keep.len() && keep[i] && keep[j] + }); + for a in &mut self.arcs { + a.n0 = remap[a.n0 as usize]; + a.n1 = remap[a.n1 as usize]; + } + + removed + } + + /// removes selected segments and returns the count. + pub fn delete_selected_segments(&mut self) -> usize { + let before = self.segments.len(); + self.segments.retain(|s| !s.selected); + before - self.segments.len() + } + + /// removes selected arc segments and returns the count. + pub fn delete_selected_arcs(&mut self) -> usize { + let before = self.arcs.len(); + self.arcs.retain(|a| !a.selected); + before - self.arcs.len() + } + + /// removes selected block labels and returns the count. + pub fn delete_selected_block_labels(&mut self) -> usize { + let before = self.block_labels.len(); + self.block_labels.retain(|b| !b.selected); + before - self.block_labels.len() + } + + /// returns the index of the closest node to (x, y), or None when the node list is empty. + pub fn closest_node(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, n) in self.nodes.iter().enumerate() { + let d = (n.x - x).hypot(n.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the closest block label to (x, y), or None when none exist. + pub fn closest_block_label(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, b) in self.block_labels.iter().enumerate() { + let d = (b.x - x).hypot(b.y - y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// returns the index of the segment whose nearest point is closest to (x, y). + pub fn closest_segment(&self, x: f64, y: f64) -> Option { + let mut best: Option<(usize, f64)> = None; + for (i, s) in self.segments.iter().enumerate() { + let (Some(p0), Some(p1)) = (self.nodes.get(s.n0 as usize), self.nodes.get(s.n1 as usize)) else { continue }; + let d = point_to_segment_distance(x, y, p0.x, p0.y, p1.x, p1.y); + match best { + None => best = Some((i, d)), + Some((_, bd)) if d < bd => best = Some((i, d)), + _ => {} + } + } + best.map(|(i, _)| i) + } + + /// clears the selection flag on every geometric entity in the doc. + pub fn clear_selection(&mut self) { + for n in &mut self.nodes { n.selected = false; } + for s in &mut self.segments { s.selected = false; } + for a in &mut self.arcs { a.selected = false; } + for b in &mut self.block_labels { b.selected = false; } + } +} + +/// auto-tolerance derived from the bounding box of a node slice. +fn nodes_bbox_tolerance(nodes: &[Node]) -> f64 { + if nodes.len() < 2 { + return 1.0e-8; + } + let mut xmin = f64::INFINITY; + let mut xmax = f64::NEG_INFINITY; + let mut ymin = f64::INFINITY; + let mut ymax = f64::NEG_INFINITY; + for n in nodes { + xmin = xmin.min(n.x); xmax = xmax.max(n.x); + ymin = ymin.min(n.y); ymax = ymax.max(n.y); + } + let dx = xmax - xmin; + let dy = ymax - ymin; + (dx * dx + dy * dy).sqrt() * BBOX_TOLERANCE_FRAC +} + +/// Euclidean distance from (px, py) to the segment between (ax, ay) and (bx, by). +fn point_to_segment_distance(px: f64, py: f64, ax: f64, ay: f64, bx: f64, by: f64) -> f64 { + let dx = bx - ax; + let dy = by - ay; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return (px - ax).hypot(py - ay); + } + let t = (((px - ax) * dx + (py - ay) * dy) / len2).clamp(0.0, 1.0); + let cx = ax + t * dx; + let cy = ay + t * dy; + (px - cx).hypot(py - cy) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn doc_with_corners() -> FemmDoc { + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0,-1.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d + } + + #[test] + fn duplicate_segment_rejected() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert!(!d.add_segment(0, 1)); + assert!(!d.add_segment(1, 0)); + assert_eq!(d.segments.len(), 1); + } + + #[test] + fn segment_passing_through_existing_node_splits() { + // three colinear nodes (-1,0), (1,0), (0,0). adding the outer-to-outer segment + // splits at the midpoint node, producing two pieces meeting there. + let mut d = FemmDoc::default(); + d.add_node(-1.0, 0.0, 0.0); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 0.0, 0.0); + assert!(d.add_segment(0, 1)); + assert_eq!(d.segments.len(), 2); + let touches_mid = d.segments.iter().filter(|s| s.n0 == 2 || s.n1 == 2).count(); + assert_eq!(touches_mid, 2); + } + + #[test] + fn enforce_pslg_splits_first_segment_at_late_intersection() { + // incremental edit splits only the second of two crossing segments. enforce_pslg + // back-splits the first segment through the intersection node. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + assert_eq!(d.segments.len(), 3, "pre-enforce: one whole, one split"); + + d.enforce_pslg(); + + assert_eq!(d.nodes.len(), 5, "no node added or merged by enforce"); + assert_eq!(d.segments.len(), 4, "both segments split at the origin"); + let touches_origin = d.segments.iter().filter(|s| s.n0 == 4 || s.n1 == 4).count(); + assert_eq!(touches_origin, 4); + } + + #[test] + fn enforce_pslg_preserves_segment_boundary_marker() { + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + assert!(d.add_segment_with_marker(0, 1, "outer")); + + d.enforce_pslg(); + + assert_eq!(d.segments.len(), 1); + assert_eq!(d.segments[0].boundary_marker, "outer"); + } + + #[test] + fn enforce_pslg_preserves_arc_metadata() { + // arc carrying non-default marker, side length, group, and direction. + // enforce_pslg round-trips every field intact. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let template = ArcSegment { + n0: 0, + n1: 1, + arc_length: 90.0, + max_side_length: 5.0, + boundary_marker: String::from("outer"), + in_conductor: String::new(), + hidden: true, + in_group: 7, + normal_direction: false, + selected: false, + }; + assert!(d.add_arc_segment_with_template(0, 1, 90.0, &template)); + + d.enforce_pslg(); + + assert_eq!(d.arcs.len(), 1); + let a = &d.arcs[0]; + assert_eq!(a.boundary_marker, "outer"); + assert!((a.max_side_length - 5.0).abs() < 1e-12); + assert!(a.hidden); + assert_eq!(a.in_group, 7); + assert!(!a.normal_direction); + } + + #[test] + fn duplicate_arc_rejected() { + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(!d.add_arc_segment(0, 1, 90.0)); + // reversed endpoints define a distinct arc. + assert!(d.add_arc_segment(1, 0, 90.0)); + assert_eq!(d.arcs.len(), 2); + } + + #[test] + fn add_segment_crossing_unit_quarter_arc_splits_segment() { + // unit quarter arc from (1,0) to (0,1); horizontal line at y=0.5 crosses the arc once. + // expect a fresh intersection node near (0.866, 0.5) and the segment split into two. + let mut d = FemmDoc::default(); + d.add_node( 1.0, 0.0, 0.0); + d.add_node( 0.0, 1.0, 0.0); + d.add_node(-2.0, 0.5, 0.0); + d.add_node( 2.0, 0.5, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5); + let n4 = &d.nodes[4]; + assert!((n4.x - 0.75_f64.sqrt()).abs() < 1e-6); + assert!((n4.y - 0.5).abs() < 1e-9); + assert_eq!(d.segments.len(), 2); + assert_eq!(d.arcs.len(), 1); + } + + #[test] + fn add_arc_passing_through_existing_node_splits_into_two_arcs() { + // pre-existing node at (cos45, sin45) lies on the unit quarter arc from (1,0) to (0,1). + // adding the arc detects the on-arc node and emits two sub-arcs of 45 degrees each. + let mut d = FemmDoc::default(); + d.add_node(1.0, 0.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + let mid_x = std::f64::consts::FRAC_1_SQRT_2; + d.add_node(mid_x, mid_x, 0.0); + assert!(d.add_arc_segment(0, 1, 90.0)); + + assert_eq!(d.arcs.len(), 2); + let touches_mid = d.arcs.iter().filter(|a| a.n0 == 2 || a.n1 == 2).count(); + assert_eq!(touches_mid, 2); + // split sweeps sum to the parent 90 degrees, within floating slop. + let total: f64 = d.arcs.iter().map(|a| a.arc_length).sum(); + assert!((total - 90.0).abs() < 1e-6, "sum of split sweeps = {total}"); + } + + #[test] + fn crossing_second_segment_splits_at_intersection() { + // two perpendicular segments through the origin. horizontal first, vertical second. + // the vertical add inserts an origin node and self-splits; the horizontal stays whole. + let mut d = doc_with_corners(); + assert!(d.add_segment(0, 1)); + assert!(d.add_segment(2, 3)); + + assert_eq!(d.nodes.len(), 5, "intersection node added"); + let new_idx = 4; + assert!((d.nodes[new_idx].x).abs() < 1e-9); + assert!((d.nodes[new_idx].y).abs() < 1e-9); + + // expected pieces: 0->1 (horizontal whole), 2->4 and 4->3 (vertical halves). + assert_eq!(d.segments.len(), 3); + let touches_new = d.segments.iter().filter(|s| s.n0 == new_idx as i32 || s.n1 == new_idx as i32).count(); + assert_eq!(touches_new, 2); + } +} diff --git a/crates/femm-doc-heat/src/geom_math.rs b/crates/femm-doc-heat/src/geom_math.rs new file mode 100644 index 0000000..37eeb54 --- /dev/null +++ b/crates/femm-doc-heat/src/geom_math.rs @@ -0,0 +1,322 @@ +//! 2D intersection math and distance queries over (x, y) pairs and degree-sweep arcs. + +use num_complex::Complex64; + +/// closeness factor for common-endpoint and small-interval rejection. +const EPS_FRAC: f64 = 1.0e-8; + +/// recovers the center and radius of the circle through an arc's endpoints with the given sweep. +pub fn circle_from_arc(p0: (f64, f64), p1: (f64, f64), sweep_deg: f64) -> (f64, f64, f64) { + let a0 = Complex64::new(p0.0, p0.1); + let a1 = Complex64::new(p1.0, p1.1); + let chord = (a1 - a0).norm(); + let t = (a1 - a0) / chord; + let tta = sweep_deg.to_radians(); + let r = chord / (2.0 * (tta / 2.0).sin()); + let c = a0 + (Complex64::new(chord / 2.0, (r * r - chord * chord / 4.0).sqrt())) * t; + (c.re, c.im, r.abs()) +} + +/// shortest distance from a point to a segment. +pub fn shortest_distance_from_segment(p: (f64, f64), a: (f64, f64), b: (f64, f64)) -> f64 { + let dx = b.0 - a.0; + let dy = b.1 - a.1; + let len2 = dx * dx + dy * dy; + if len2 < 1e-18 { + return ((p.0 - a.0).powi(2) + (p.1 - a.1).powi(2)).sqrt(); + } + let t = (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / len2).clamp(0.0, 1.0); + let cx = a.0 + t * dx; + let cy = a.1 + t * dy; + ((p.0 - cx).powi(2) + (p.1 - cy).powi(2)).sqrt() +} + +/// shortest distance from a point to a sweep-defined arc. +pub fn shortest_distance_from_arc( + p: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> f64 { + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let pc = Complex64::new(p.0 - cx, p.1 - cy); + let d = pc.norm(); + if d == 0.0 { + return r; + } + let t = pc / d; + let foot = pc - Complex64::new(r * t.re, r * t.im); + let l = foot.norm(); + let a0_dir = Complex64::new(a0.0 - cx, a0.1 - cy); + let z = ((t / a0_dir).arg() * 180.0 / std::f64::consts::PI + 360.0) % 360.0; + let sweep_abs = sweep_deg.abs(); + if z > 0.0 && z < sweep_abs { + return l; + } + let e0 = ((p.0 - a0.0).powi(2) + (p.1 - a0.1).powi(2)).sqrt(); + let e1 = ((p.0 - a1.0).powi(2) + (p.1 - a1.1).powi(2)).sqrt(); + e0.min(e1) +} + +/// intersection of the open segments p0->p1 and q0->q1, or None when none exists. +pub fn line_line_intersection( + p0: (f64, f64), + p1: (f64, f64), + q0: (f64, f64), + q1: (f64, f64), +) -> Option<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let q0c = Complex64::new(q0.0, q0.1); + let q1c = Complex64::new(q1.0, q1.1); + + // shared endpoint -> no other intersection + if (p0c - q0c).norm() < f64::EPSILON + || (p0c - q1c).norm() < f64::EPSILON + || (p1c - q0c).norm() < f64::EPSILON + || (p1c - q1c).norm() < f64::EPSILON + { + return None; + } + + let ee = ((p1c - p0c).norm()).min((q1c - q0c).norm()) * EPS_FRAC; + + let denom = p1c - p0c; + if denom.norm() < f64::EPSILON { + return None; + } + let r0 = (q0c - p0c) / denom; + let r1 = (q1c - p0c) / denom; + + if r0.re <= 0.0 && r1.re <= 0.0 { return None; } + if r0.re >= 1.0 && r1.re >= 1.0 { return None; } + if r0.im <= 0.0 && r1.im <= 0.0 { return None; } + if r0.im >= 0.0 && r1.im >= 0.0 { return None; } + + let denom_im = r0.im - r1.im; + if denom_im.abs() < f64::EPSILON { return None; } + let z = r0.im / denom_im; + let x = ((1.0 - z) * r0 + z * r1).re; + if x < ee || x > 1.0 - ee { + return None; + } + + let hit = Complex64::new(p0.0, p0.1) * (1.0 - z) + Complex64::new(q0.0, q0.1) * 0.0; + let _ = hit; + let result = Complex64::new(q0.0, q0.1) * (1.0 - z) + Complex64::new(q1.0, q1.1) * z; + Some((result.re, result.im)) +} + +/// intersection points of an open line segment with an open arc, returning 0, 1, or 2 hits. +pub fn line_arc_intersection( + p0: (f64, f64), + p1: (f64, f64), + a0: (f64, f64), + a1: (f64, f64), + sweep_deg: f64, +) -> Vec<(f64, f64)> { + let p0c = Complex64::new(p0.0, p0.1); + let p1c = Complex64::new(p1.0, p1.1); + let a0c = Complex64::new(a0.0, a0.1); + let (cx, cy, r) = circle_from_arc(a0, a1, sweep_deg); + let c = Complex64::new(cx, cy); + + let d = (p1c - p0c).norm(); + if d < f64::EPSILON { + return Vec::new(); + } + let t = (p1c - p0c) / d; + let v = (c - p0c) / t; + + if v.im.abs() > r { + return Vec::new(); + } + let l = (r * r - v.im * v.im).sqrt(); + let tta = sweep_deg.to_radians().abs(); + + let mut out = Vec::with_capacity(2); + + if (l / r) < 1.0e-5 { + let hit = p0c + Complex64::new(v.re, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + return out; + } + + for sign in [1.0, -1.0] { + let hit = p0c + Complex64::new(v.re + sign * l, 0.0) * t; + let r_param = ((hit - p0c) / t).re; + let z = ((hit - c) / (a0c - c)).arg(); + if r_param > 0.0 && r_param < d && z > 0.0 && z < tta { + out.push((hit.re, hit.im)); + } + } + out +} + +/// intersection points of two open arcs, returning 0, 1, or 2 hits. +pub fn arc_arc_intersection( + a0_p0: (f64, f64), + a0_p1: (f64, f64), + a0_sweep_deg: f64, + a1_p0: (f64, f64), + a1_p1: (f64, f64), + a1_sweep_deg: f64, +) -> Vec<(f64, f64)> { + let (c0x, c0y, r0) = circle_from_arc(a0_p0, a0_p1, a0_sweep_deg); + let (c1x, c1y, r1) = circle_from_arc(a1_p0, a1_p1, a1_sweep_deg); + let c0 = Complex64::new(c0x, c0y); + let c1 = Complex64::new(c1x, c1y); + let d = (c1 - c0).norm(); + + if d > r0 + r1 || d < 1.0e-8 { + return Vec::new(); + } + + let l = ((r0 + r1 - d) * (d + r0 - r1) * (d - r0 + r1) * (d + r0 + r1)).sqrt() / (2.0 * d); + let c = 1.0 + (r0 / d) * (r0 / d) - (r1 / d) * (r1 / d); + let t = (c1 - c0) / d; + let tta0 = a0_sweep_deg.to_radians().abs(); + let tta1 = a1_sweep_deg.to_radians().abs(); + let a0c = Complex64::new(a0_p0.0, a0_p0.1); + let a1c = Complex64::new(a1_p0.0, a1_p0.1); + + let mut out = Vec::with_capacity(2); + + let first = c0 + Complex64::new(c * d / 2.0, l) * t; + let z0 = ((first - c0) / (a0c - c0)).arg(); + let z1 = ((first - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((first.re, first.im)); + } + + // tangent-touch case: only one intersection + if (d - r0 + r1).abs() / (r0 + r1) < 1.0e-5 { + return out; + } + + let second = c0 + Complex64::new(c * d / 2.0, -l) * t; + let z0 = ((second - c0) / (a0c - c0)).arg(); + let z1 = ((second - c1) / (a1c - c1)).arg(); + if z0 > 0.0 && z0 < tta0 && z1 > 0.0 && z1 < tta1 { + out.push((second.re, second.im)); + } + + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: (f64, f64), b: (f64, f64), tol: f64) -> bool { + ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt() < tol + } + + #[test] + fn line_line_crosses_at_origin() { + // (-1,-1)->(1,1) crosses (-1,1)->(1,-1) at origin. + let hit = line_line_intersection((-1.0, -1.0), (1.0, 1.0), (-1.0, 1.0), (1.0, -1.0)); + assert!(hit.is_some()); + assert!(close(hit.unwrap(), (0.0, 0.0), 1e-9)); + } + + #[test] + fn line_line_parallel_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn line_line_t_junction_at_endpoint_rejects() { + // the second segment's q0 sits ON the first segment's interior + // but the prospective endpoint check should reject (intersection at q0). + let hit = line_line_intersection((0.0, 0.0), (2.0, 0.0), (1.0, 0.0), (1.0, 1.0)); + // shape: one endpoint sits on the other segment, no clean crossing. + // GetIntersection returns FALSE when an intersection lies within ee of an endpoint of the + // prospective line — here the meeting point IS the prospective q0, so no split is wanted. + assert!(hit.is_none()); + } + + #[test] + fn line_line_shared_endpoint_no_intersection() { + let hit = line_line_intersection((0.0, 0.0), (1.0, 0.0), (1.0, 0.0), (2.0, 1.0)); + assert!(hit.is_none()); + } + + #[test] + fn circle_from_quarter_arc_unit() { + // a quarter circle (90 deg sweep) from (1, 0) to (0, 1) sits on the unit circle centered at origin. + let (cx, cy, r) = circle_from_arc((1.0, 0.0), (0.0, 1.0), 90.0); + assert!((cx - 0.0).abs() < 1e-9); + assert!((cy - 0.0).abs() < 1e-9); + assert!((r - 1.0).abs() < 1e-9); + } + + #[test] + fn line_through_quarter_arc_crosses_once() { + // diagonal line from (-2, 0.5) to (2, 0.5) crosses the quarter circle (1,0)->(0,1) at one point. + let hits = line_arc_intersection((-2.0, 0.5), (2.0, 0.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert_eq!(hits.len(), 1); + let h = hits[0]; + assert!((h.0 * h.0 + h.1 * h.1 - 1.0).abs() < 1e-9, "point should lie on unit circle: {h:?}"); + assert!((h.1 - 0.5).abs() < 1e-9); + assert!(h.0 > 0.0); + } + + #[test] + fn line_outside_arc_does_not_cross() { + // line clearly outside the arc swept region. + let hits = line_arc_intersection((-2.0, 2.5), (2.0, 2.5), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!(hits.is_empty()); + } + + #[test] + fn arc_arc_two_intersections() { + // two unit circles offset by 1.0 along x, both swept 360 degrees... but we use partial arcs + // shaped to clip both crossings in their swept range. + // circle A: center (0,0), arcs covering full top half (180 deg from (1,0) to (-1,0)). + // circle B: center (1,0), arcs covering full top half (180 deg from (2,0) to (0,0)). + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (2.0, 0.0), (0.0, 0.0), 180.0, + ); + assert_eq!(hits.len(), 1, "two top-half arcs of overlapping unit circles meet at one upper crossing: {hits:?}"); + let h = hits[0]; + assert!((h.0 - 0.5).abs() < 1e-9, "x = 0.5 by symmetry, got {h:?}"); + assert!(h.1 > 0.0, "intersection should be above the x-axis: {h:?}"); + } + + #[test] + fn arc_arc_disjoint() { + let hits = arc_arc_intersection( + (1.0, 0.0), (-1.0, 0.0), 180.0, + (11.0, 0.0), (9.0, 0.0), 180.0, + ); + assert!(hits.is_empty()); + } + + #[test] + fn shortest_distance_from_segment_perpendicular() { + // point (0, 1) to segment (-1, 0)->(1, 0) is distance 1. + let d = shortest_distance_from_segment((0.0, 1.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_segment_clamps_to_endpoint() { + // point (-2, 0) to segment (-1, 0)->(1, 0) is distance 1 (clamps to left endpoint). + let d = shortest_distance_from_segment((-2.0, 0.0), (-1.0, 0.0), (1.0, 0.0)); + assert!((d - 1.0).abs() < 1e-9); + } + + #[test] + fn shortest_distance_from_arc_radial() { + // point at origin distance to unit-circle quarter arc is 1. + let d = shortest_distance_from_arc((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), 90.0); + assert!((d - 1.0).abs() < 1e-9); + } +} diff --git a/crates/femm-doc-heat/src/lib.rs b/crates/femm-doc-heat/src/lib.rs index 2fa5ab2..cfbead4 100644 --- a/crates/femm-doc-heat/src/lib.rs +++ b/crates/femm-doc-heat/src/lib.rs @@ -1,9 +1,11 @@ //! heat-flow pre-processor document model: geometry, properties, parser, writer. pub mod geom; +pub mod geom_math; pub mod props; pub mod parser; pub mod writer; +pub mod edit; use num_complex::Complex64; diff --git a/crates/femm-doc-mag/src/edit.rs b/crates/femm-doc-mag/src/edit.rs index 0fd5b3f..2b92fd1 100644 --- a/crates/femm-doc-mag/src/edit.rs +++ b/crates/femm-doc-mag/src/edit.rs @@ -553,8 +553,8 @@ mod tests { #[test] fn add_segment_crossing_unit_quarter_arc_splits_segment() { - // unit quarter arc from (1,0) to (0,1); horizontal line at y=0.5 crosses it once. - // expect a new intersection node near (0.866, 0.5) and the segment split into two. + // unit quarter arc from (1,0) to (0,1); horizontal line at y=0.5 crosses the arc once. + // expect a fresh intersection node near (0.866, 0.5) and the segment split into two. let mut d = FemmDoc::default(); d.add_node( 1.0, 0.0, 0.0); d.add_node( 0.0, 1.0, 0.0); @@ -574,7 +574,7 @@ mod tests { #[test] fn add_arc_passing_through_existing_node_splits_into_two_arcs() { // pre-existing node at (cos45, sin45) lies on the unit quarter arc from (1,0) to (0,1). - // adding the arc should detect the on-arc node and emit two sub-arcs of 45 degrees each. + // adding the arc detects the on-arc node and emits two sub-arcs of 45 degrees each. let mut d = FemmDoc::default(); d.add_node(1.0, 0.0, 0.0); d.add_node(0.0, 1.0, 0.0); @@ -592,9 +592,8 @@ mod tests { #[test] fn crossing_second_segment_splits_at_intersection() { - // horizontal (-1,0)->(1,0) added first stays whole. - // the vertical (0,-1)->(0,1) is added second; it intersects at origin, - // so a new node lands at (0,0) and the vertical splits into two pieces. + // two perpendicular segments through the origin. horizontal first, vertical second. + // the vertical add inserts an origin node and self-splits; the horizontal stays whole. let mut d = doc_with_corners(); assert!(d.add_segment(0, 1)); assert!(d.add_segment(2, 3)); @@ -604,8 +603,7 @@ mod tests { assert!((d.nodes[new_idx].x).abs() < 1e-9); assert!((d.nodes[new_idx].y).abs() < 1e-9); - // exactly one piece touches node 0 (-1,0) and node 4 (0,0): the original horizontal stays whole here. - // expected pieces: 0->1 (horizontal whole), 2->4, 4->3 (vertical halves). + // expected pieces: 0->1 (horizontal whole), 2->4 and 4->3 (vertical halves). assert_eq!(d.segments.len(), 3); let touches_new = d.segments.iter().filter(|s| s.n0 == new_idx as i32 || s.n1 == new_idx as i32).count(); assert_eq!(touches_new, 2); diff --git a/crates/femm-doc-mag/src/lib.rs b/crates/femm-doc-mag/src/lib.rs index 9143cd3..7c931e2 100644 --- a/crates/femm-doc-mag/src/lib.rs +++ b/crates/femm-doc-mag/src/lib.rs @@ -6,6 +6,7 @@ pub mod props; pub mod parser; pub mod writer; pub mod edit; +pub mod poly; use num_complex::Complex64; diff --git a/crates/femm-doc-mag/src/poly.rs b/crates/femm-doc-mag/src/poly.rs new file mode 100644 index 0000000..925e86e --- /dev/null +++ b/crates/femm-doc-mag/src/poly.rs @@ -0,0 +1,166 @@ +//! .poly emitter consumed by the Triangle mesher. + +use crate::FemmDoc; +use std::fmt::Write; + +impl FemmDoc { + /// renders the doc geometry to the Triangle .poly text format, returning the file contents. + pub fn write_poly(&self) -> String { + let mut out = String::new(); + let bbox_diag = self.bbox_diagonal(); + let default_area = if bbox_diag > 0.0 { bbox_diag } else { -1.0 }; + + writeln!(out, "{}\t2\t0\t1", self.nodes.len()).unwrap(); + for (i, n) in self.nodes.iter().enumerate() { + let marker = point_marker_index(self, &n.boundary_marker); + writeln!(out, "{i}\t{:.17e}\t{:.17e}\t{marker}", n.x, n.y).unwrap(); + } + + writeln!(out, "{}\t1", self.segments.len()).unwrap(); + for (i, s) in self.segments.iter().enumerate() { + let marker = -boundary_marker_index(self, &s.boundary_marker); + writeln!(out, "{i}\t{}\t{}\t{marker}", s.n0, s.n1).unwrap(); + } + + let holes: Vec<&_> = self.block_labels.iter().filter(|b| b.block_type == "").collect(); + writeln!(out, "{}", holes.len()).unwrap(); + for (i, h) in holes.iter().enumerate() { + writeln!(out, "{i}\t{:.17e}\t{:.17e}", h.x, h.y).unwrap(); + } + + let regions: Vec<&_> = self.block_labels.iter().filter(|b| b.block_type != "").collect(); + writeln!(out, "{}", regions.len()).unwrap(); + for (i, r) in regions.iter().enumerate() { + let attr = (i as i32) + 1; + let max_area = if r.max_area > 0.0 && r.max_area < default_area { + r.max_area + } else { + default_area + }; + writeln!(out, "{i}\t{:.17e}\t{:.17e}\t{attr}\t{:.17e}", r.x, r.y, max_area).unwrap(); + } + + out + } + + /// writes the .poly text to disk at the given path. + pub fn save_poly(&self, path: impl AsRef) -> std::io::Result<()> { + std::fs::write(path, self.write_poly()) + } + + fn bbox_diagonal(&self) -> f64 { + if self.nodes.len() < 2 { return 0.0; } + let mut xmin = f64::INFINITY; + let mut xmax = f64::NEG_INFINITY; + let mut ymin = f64::INFINITY; + let mut ymax = f64::NEG_INFINITY; + for n in &self.nodes { + xmin = xmin.min(n.x); xmax = xmax.max(n.x); + ymin = ymin.min(n.y); ymax = ymax.max(n.y); + } + let dx = xmax - xmin; + let dy = ymax - ymin; + (dx * dx + dy * dy).sqrt() + } +} + +/// resolves a node's boundary-marker name to Triangle's marker integer (0 for unmarked). +fn point_marker_index(doc: &FemmDoc, name: &str) -> i32 { + if name.is_empty() { return 0; } + for (i, p) in doc.points.iter().enumerate() { + if p.name == name { return (i as i32) + 2; } + } + 0 +} + +/// resolves a segment's boundary-marker name to Triangle's marker integer, returned positive. +/// the caller negates for Triangle's segment-marker convention. +fn boundary_marker_index(doc: &FemmDoc, name: &str) -> i32 { + if name.is_empty() { return 0; } + for (i, b) in doc.boundaries.iter().enumerate() { + if b.name == name { return (i as i32) + 2; } + } + 0 +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{BoundaryProp, PointProp, BlockLabel}; + + #[test] + fn write_poly_header_lines_count_correctly() { + // four-corner square with one block label inside; verify the section headers and counts. + let mut d = FemmDoc::default(); + d.points.push(PointProp { name: "A=0".into(), ..PointProp::default() }); + d.boundaries.push(BoundaryProp { name: "outer".into(), ..BoundaryProp::default() }); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + d.add_node(1.0, 1.0, 0.0); + d.add_node(0.0, 1.0, 0.0); + d.nodes[0].boundary_marker = "A=0".into(); + d.add_segment_with_marker(0, 1, "outer"); + d.add_segment_with_marker(1, 2, "outer"); + d.add_segment_with_marker(2, 3, "outer"); + d.add_segment_with_marker(3, 0, "outer"); + d.block_labels.push(BlockLabel { + x: 0.5, y: 0.5, + max_area: 0.0, + block_type: "Air".into(), + ..BlockLabel::default() + }); + + let poly = d.write_poly(); + let lines: Vec<&str> = poly.lines().collect(); + + assert!(lines[0].starts_with("4\t2\t0\t1"), "node header: {}", lines[0]); + assert!(lines[5].starts_with("4\t1"), "segment header: {}", lines[5]); + assert_eq!(lines[10], "0", "no holes section: {}", lines[10]); + assert_eq!(lines[11], "1", "one region: {}", lines[11]); + } + + #[test] + fn write_poly_marks_nodes_and_segments_by_property_index() { + // first point property listed -> marker 2; first boundary property listed -> -2 on segment. + let mut d = FemmDoc::default(); + d.points.push(PointProp { name: "A=0".into(), ..PointProp::default() }); + d.boundaries.push(BoundaryProp { name: "outer".into(), ..BoundaryProp::default() }); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + d.nodes[0].boundary_marker = "A=0".into(); + d.add_segment_with_marker(0, 1, "outer"); + + let poly = d.write_poly(); + let lines: Vec<&str> = poly.lines().collect(); + + assert!(lines[1].ends_with("\t2"), "node 0 marker = 2: {}", lines[1]); + assert!(lines[2].ends_with("\t0"), "node 1 marker = 0: {}", lines[2]); + assert!(lines[4].ends_with("\t-2"), "segment 0 marker = -2: {}", lines[4]); + } + + #[test] + fn write_poly_holes_partitioned_from_regions() { + // one "" label counts as a hole; everything else is a region. + let mut d = FemmDoc::default(); + d.add_node(0.0, 0.0, 0.0); + d.add_node(1.0, 0.0, 0.0); + d.block_labels.push(BlockLabel { + x: 0.5, y: 0.5, + block_type: "".into(), + ..BlockLabel::default() + }); + d.block_labels.push(BlockLabel { + x: 0.2, y: 0.2, + block_type: "Air".into(), + ..BlockLabel::default() + }); + + let poly = d.write_poly(); + let lines: Vec<&str> = poly.lines().collect(); + + // node header (0), 2 node rows (1,2), segment header "0\t1" (3), + // hole count (4), 1 hole row (5), region count (6), 1 region row (7). + assert_eq!(lines[4], "1", "one hole: {}", lines[4]); + assert_eq!(lines[6], "1", "one region: {}", lines[6]); + } +} diff --git a/scripts/macos/build.sh b/scripts/macos/build.sh new file mode 100755 index 0000000..da5cde7 --- /dev/null +++ b/scripts/macos/build.sh @@ -0,0 +1,86 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT="$(cd "$(dirname "$0")/../.." && pwd)" +cd "$ROOT" + +case "$(uname -s)" in + Darwin) ;; + *) echo "wrong platform: $(uname -s) - use cargo xtask build" >&2; exit 1 ;; +esac + +export MACOSX_DEPLOYMENT_TARGET=11.0 + +BUILD="$ROOT/build" +APP="$BUILD/bin/femm.app" +CONTENTS="$APP/Contents" +MACOS="$CONTENTS/MacOS" +RESOURCES="$CONTENTS/Resources" +ICONSET="$BUILD/AppIcon.iconset" +ICNS="$RESOURCES/AppIcon.icns" +SVG="$ROOT/assets/femm.svg" + +echo "Building Rust workspace (release)..." +cargo build --release -p femm-app + +BIN="$ROOT/target/release/femm" +if [ ! -f "$BIN" ]; then + echo "ERROR: femm binary not found at $BIN" >&2 + exit 1 +fi + +rm -rf "$APP" +mkdir -p "$MACOS" "$RESOURCES" +cp "$BIN" "$MACOS/femm" + +if [ -f "$SVG" ]; then + if ! command -v rsvg-convert >/dev/null; then + echo "ERROR: rsvg-convert missing - brew install librsvg" >&2 + exit 1 + fi + if ! command -v iconutil >/dev/null; then + echo "ERROR: iconutil missing - install Xcode command line tools" >&2 + exit 1 + fi + rm -rf "$ICONSET" + mkdir -p "$ICONSET" + for size in 16 32 64 128 256 512 1024; do + rsvg-convert --width="$size" --height="$size" "$SVG" -o "$ICONSET/icon_${size}.png" + done + cp "$ICONSET/icon_16.png" "$ICONSET/icon_16x16.png" + cp "$ICONSET/icon_32.png" "$ICONSET/icon_16x16@2x.png" + cp "$ICONSET/icon_32.png" "$ICONSET/icon_32x32.png" + cp "$ICONSET/icon_64.png" "$ICONSET/icon_32x32@2x.png" + cp "$ICONSET/icon_128.png" "$ICONSET/icon_128x128.png" + cp "$ICONSET/icon_256.png" "$ICONSET/icon_128x128@2x.png" + cp "$ICONSET/icon_256.png" "$ICONSET/icon_256x256.png" + cp "$ICONSET/icon_512.png" "$ICONSET/icon_256x256@2x.png" + cp "$ICONSET/icon_512.png" "$ICONSET/icon_512x512.png" + cp "$ICONSET/icon_1024.png" "$ICONSET/icon_512x512@2x.png" + iconutil -c icns "$ICONSET" -o "$ICNS" +else + echo "WARNING: $SVG not found - app bundle will lack an icon" +fi + +cat > "$CONTENTS/Info.plist" <<'PLIST' + + + + + CFBundleDevelopmentRegion en + CFBundleExecutable femm + CFBundleIconFile AppIcon + CFBundleIdentifier org.else-if.femm + CFBundleInfoDictionaryVersion 6.0 + CFBundleName femm + CFBundleDisplayName FEMM + CFBundlePackageType APPL + CFBundleShortVersionString 0.0.1 + CFBundleVersion 0.0.1 + LSMinimumSystemVersion 11.0 + NSHighResolutionCapable + + +PLIST + +echo "Built: $APP" diff --git a/scripts/macos/debug.sh b/scripts/macos/debug.sh new file mode 100755 index 0000000..b7d5f26 --- /dev/null +++ b/scripts/macos/debug.sh @@ -0,0 +1,13 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT="$(cd "$(dirname "$0")/../.." && pwd)" +cd "$ROOT" + +case "$(uname -s)" in + Darwin) ;; + *) echo "wrong platform: $(uname -s) - use cargo xtask debug" >&2; exit 1 ;; +esac + +export RUST_BACKTRACE=1 +exec cargo run -p femm-app -- "$@" diff --git a/scripts/macos/ffi.sh b/scripts/macos/ffi.sh new file mode 100755 index 0000000..12be94a --- /dev/null +++ b/scripts/macos/ffi.sh @@ -0,0 +1,5 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT="$(cd "$(dirname "$0")/../.." && pwd)" +exec bash "$ROOT/scripts/macos/build_ffi.sh" "$@" diff --git a/scripts/macos/install.sh b/scripts/macos/install.sh new file mode 100755 index 0000000..1c4e5e9 --- /dev/null +++ b/scripts/macos/install.sh @@ -0,0 +1,23 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT="$(cd "$(dirname "$0")/../.." && pwd)" +cd "$ROOT" + +case "$(uname -s)" in + Darwin) ;; + *) echo "wrong platform: $(uname -s) - use cargo xtask install" >&2; exit 1 ;; +esac + +DEST="/Applications/femm.app" + +bash "$ROOT/scripts/macos/build.sh" + +pkill -f "femm.app/Contents/MacOS/femm" 2>/dev/null || true +sleep 0.5 + +echo "Installing to $DEST..." +rm -rf "$DEST" +cp -R "$ROOT/build/bin/femm.app" "$DEST" + +echo "Installed: $DEST" diff --git a/xtask/Cargo.toml b/xtask/Cargo.toml new file mode 100644 index 0000000..f706b5b --- /dev/null +++ b/xtask/Cargo.toml @@ -0,0 +1,9 @@ +[package] +name = "xtask" +version = "0.0.0" +edition = "2021" +publish = false + +[[bin]] +name = "xtask" +path = "src/main.rs" diff --git a/xtask/src/main.rs b/xtask/src/main.rs new file mode 100644 index 0000000..3833458 --- /dev/null +++ b/xtask/src/main.rs @@ -0,0 +1,119 @@ +use std::env; +use std::path::PathBuf; +use std::process::{Command, ExitCode}; + +const KNOWN_PLATFORMS: &[&str] = &["macos", "windows", "linux"]; + +/// dispatches a cargo xtask sub-command to the matching platform script under scripts/. +fn main() -> ExitCode { + let args: Vec = env::args().skip(1).collect(); + let cmd = args.first().map(String::as_str).unwrap_or(""); + + if cmd.is_empty() || cmd == "help" || cmd == "--help" || cmd == "-h" { + print_help(); + return ExitCode::from(2); + } + + let extra_args: Vec<&String> = args.iter().skip(1).collect(); + let (action, platform) = parse(cmd); + + let repo_root = PathBuf::from(env!("CARGO_MANIFEST_DIR")) + .parent() + .expect("xtask manifest must have a parent") + .to_path_buf(); + + let (script, runner) = match platform.as_str() { + "windows" => ( + repo_root.join(format!("scripts/windows/{action}.ps1")), + vec![ + "powershell", + "-NoProfile", + "-ExecutionPolicy", + "Bypass", + "-File", + ], + ), + "linux" | "macos" => ( + repo_root.join(format!("scripts/{platform}/{action}.sh")), + vec!["bash"], + ), + other => { + eprintln!("unknown platform: {other}"); + return ExitCode::from(2); + } + }; + + if !script.exists() { + eprintln!("script not found: {}", script.display()); + return ExitCode::from(1); + } + + let extra_display = if extra_args.is_empty() { + String::new() + } else { + format!( + " {}", + extra_args.iter().map(|s| s.as_str()).collect::>().join(" "), + ) + }; + eprintln!("-> {} {}{}", runner.join(" "), script.display(), extra_display); + + let mut command = Command::new(runner[0]); + for arg in &runner[1..] { + command.arg(arg); + } + command.arg(&script); + for a in &extra_args { + command.arg(a.as_str()); + } + command.current_dir(&repo_root); + + match command.status() { + Ok(status) if status.success() => ExitCode::SUCCESS, + Ok(status) => ExitCode::from(status.code().unwrap_or(1) as u8), + Err(e) => { + eprintln!("failed to run {}: {e}", script.display()); + ExitCode::from(1) + } + } +} + +/// splits an action-platform command into (action, platform), defaulting to the host OS. +fn parse(cmd: &str) -> (String, String) { + if let Some(idx) = cmd.rfind('-') { + let suffix = &cmd[idx + 1..]; + if KNOWN_PLATFORMS.contains(&suffix) { + return (cmd[..idx].to_string(), suffix.to_string()); + } + } + (cmd.to_string(), current_platform().to_string()) +} + +/// resolves the host OS to a script directory name, exiting if the OS lacks scripts. +fn current_platform() -> &'static str { + match env::consts::OS { + "linux" => "linux", + "macos" => "macos", + "windows" => "windows", + other => { + eprintln!("unsupported OS: {other}"); + std::process::exit(2); + } + } +} + +/// prints command list, platform suffix syntax, and packaging hints. +fn print_help() { + eprintln!("usage: cargo xtask "); + eprintln!(); + eprintln!("commands:"); + eprintln!(" build release build for the current platform"); + eprintln!(" install release build + install (macOS: /Applications)"); + eprintln!(" debug debug build + foreground launch"); + eprintln!(" ffi rebuild the C/C++ engine archives under build/ffi/"); + eprintln!(); + eprintln!("append -macos / -windows / -linux to force a platform."); + eprintln!(" e.g. cargo xtask build-macos"); + eprintln!(); + eprintln!("trailing args pass through to the script."); +}