geop_ops_surface/
extend.rs

1//! [`ExtendSurface`]: a face standing on its own carried on past one of its
2//! edges.
3//!
4//! The face's surface is continued past the edge as the polynomial it
5//! already is (see [`NurbSurface3D::extended`]) — the natural extension: a
6//! plane stays a plane, a cylinder a cylinder, a circle carries on round,
7//! a free-form surface goes on curving as it curved. The face is rebuilt
8//! on the longer surface, bounded by its four sides.
9//!
10//! How far: `distance` along the surface, measured across the middle of
11//! the edge — the iso-curve there runs on that far, however unevenly the
12//! parametrization runs along it (see [`parameter_for_length`]). Where the
13//! parametrization runs alike all along the edge — a plane, a cylinder
14//! along its axis, any straight extrusion — the whole edge moves that far;
15//! on a free-form surface its ends move as far as the same parameter
16//! carries them.
17//!
18//! Supported: a face that is its surface's whole patch — bounded by the
19//! patch's four sides, as a boundary surface or an extruded face is — and
20//! stands alone in its sheet. Anything else is refused, by name.
21
22use geop_core_editor::{
23    edit_form::{EditForm, FormState},
24    editable::{Context, OperationGroup, OperationInfo, OperationTier},
25    editing::{Number, Unit},
26    part_operation,
27    part_operation::PartOperation,
28};
29use geop_core_geometry::{
30    contains::surface::boundary_curve, nurb_curve::NurbCurve2D, nurb_surface::NurbSurface3D,
31};
32use geop_core_math::{
33    geop_error::{GeopError, GeopResult, WithContext},
34    quadrature::{Quadrature, integrate},
35    scalars::Scalar,
36    vector::Vector2,
37    with_context,
38};
39use geop_core_program::{
40    element::ElementReference, part::StepContext, program::validate_operation_id, reader::Reader,
41};
42use geop_core_topology::{
43    Body, CoedgeGeometry, FaceId, Sense,
44    boundary::BoundaryType,
45    build::{BodySpec, CoedgeOn, CoedgeSpec, EdgeSpec, FaceSpec},
46};
47use geop_ops_model::{BodyNames, Edge, ModelOperation, NamedModel, build_model, model};
48use serde::{Deserialize, Serialize};
49
50/// Carries the face standing on its own that the edge named `edge` bounds
51/// on past that edge by `distance` (see the module docs). The face, its
52/// edges and its vertices keep their names, moved and stretched as the
53/// face grows.
54#[derive(Clone, Copy, Debug, Default, PartialEq, Serialize, Deserialize)]
55pub struct ExtendSurface;
56
57#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
58pub struct ExtendSurfaceArgs {
59    /// The edge of a face standing on its own to extend the face past.
60    pub edge: ElementReference<NamedModel, Edge>,
61    /// How far, along the surface across the middle of the edge.
62    pub distance: f64,
63}
64
65part_operation!(ExtendSurface(ExtendSurfaceArgs));
66
67impl PartOperation for ExtendSurface {
68    const KIND: &'static str = "extend_surface";
69    const INFO: OperationInfo = OperationInfo::new(
70        "Extend surface",
71        "Carry a face standing on its own on past one of its edges.",
72        OperationGroup::Surface,
73        OperationTier::Menu,
74    )
75    .experimental()
76    .icon(crate::icons::EXTEND_SURFACE);
77    type Args = ExtendSurfaceArgs;
78    type Session = ();
79
80    /// Nothing picked yet, a tenth on.
81    fn new_args(_: &Context) -> ExtendSurfaceArgs {
82        ExtendSurfaceArgs {
83            edge: ElementReference::new(model(), ""),
84            distance: 0.1,
85        }
86    }
87
88    fn run(ctx: &mut StepContext, args: &ExtendSurfaceArgs) -> GeopResult<()> {
89        let _step = ctx.step();
90        build_model::<Self>(ctx, args)
91    }
92
93    /// The edge, picked, and the distance.
94    fn form<'a>(
95        context: &'a Context<'a>,
96        args: &ExtendSurfaceArgs,
97        _: &FormState<()>,
98    ) -> EditForm<'a, ExtendSurfaceArgs> {
99        let mut f = EditForm::<ExtendSurfaceArgs>::default();
100        f.pick_one(
101            "edge",
102            "edge",
103            &context.before,
104            Some(&args.edge),
105            |edit, edge| edit.args.edge = edge,
106        );
107        let distance = Number::new("distance", args.distance, Unit::Length).range(0.0, 1.0);
108        f.number("distance", distance, |edit, d| edit.args.distance = d);
109        f
110    }
111}
112
113impl ModelOperation for ExtendSurface {
114    fn build(
115        model: &mut NamedModel,
116        _: &Reader,
117        step: &str,
118        args: &ExtendSurfaceArgs,
119    ) -> GeopResult<()> {
120        validate_operation_id(step)?;
121        if !args.distance.is_finite() {
122            return Err(GeopError::new("the distance is not a number"));
123        }
124        extend(model, args.edge.name(), Scalar::from_f64(args.distance))?;
125        Ok(())
126    }
127}
128
129/// Which side of its patch a face's coedge runs along, numbered as
130/// [`NurbSurface3D::coons`] numbers them: `v` lowest, `u` highest, `v`
131/// highest, `u` lowest — the way an outer loop runs round its patch.
132fn side_of(
133    surface: &NurbSurface3D,
134    pcurve: &geop_core_topology::Curve2,
135) -> GeopResult<Option<usize>> {
136    let ((u0, u1), (v0, v1)) = (surface.domain_u(), surface.domain_v());
137    let (a, b) = pcurve.ends()?;
138    let both = |k: usize, x: Scalar| a[k].could_be_equal(x) && b[k].could_be_equal(x);
139    // Straight along the side, not merely ending on it.
140    let straight = pcurve.degree == 1 && pcurve.control_points.len() == 2;
141    Ok(match () {
142        _ if !straight => None,
143        _ if both(1, v0) && a[0].could_be_equal(u0) && b[0].could_be_equal(u1) => Some(0),
144        _ if both(0, u1) && a[1].could_be_equal(v0) && b[1].could_be_equal(v1) => Some(1),
145        _ if both(1, v1) && a[0].could_be_equal(u1) && b[0].could_be_equal(u0) => Some(2),
146        _ if both(0, u0) && a[1].could_be_equal(v1) && b[1].could_be_equal(v0) => Some(3),
147        _ => None,
148    })
149}
150
151/// Newton steps solving for how far in the parameter an extension goes.
152const LENGTH_STEPS: usize = 30;
153
154/// How far in its `u` (`along_u`) or `v` parameter `surface` is carried on
155/// past its end (`at_end`) or start — by `carried(by)` — for the iso-curve
156/// across `(u, v)` on its border to run on `distance` past it: Newton on
157/// that length, its derivative the speed at the far end. How much to add
158/// is a free choice the distance only describes, so the steps are sharp;
159/// the length is integrated on the carried surface, as the polynomial the
160/// surface goes on as.
161fn parameter_for_length(
162    surface: &NurbSurface3D,
163    along_u: bool,
164    at_end: bool,
165    (u, v): (Scalar, Scalar),
166    distance: Scalar,
167    carried: impl Fn(Scalar) -> GeopResult<NurbSurface3D>,
168) -> GeopResult<Scalar> {
169    let speed = |s: &NurbSurface3D, t: Scalar| -> GeopResult<Scalar> {
170        let (u, v) = if along_u { (t, v) } else { (u, t) };
171        let (du, dv) = s.derivatives(u.sharpen(), v.sharpen())?;
172        Ok(if along_u { du } else { dv }.norm())
173    };
174    let from = if along_u { u } else { v };
175    // Past the end the parameter grows, before the start it falls.
176    let to = |by: Scalar| if at_end { from.add(by) } else { from.sub(by) };
177    let mut by = distance.div(speed(surface, from)?)?.sharpen();
178    let mut last = None;
179    for _ in 0..LENGTH_STEPS {
180        let longer = carried(by)?;
181        let (a, b) = if at_end {
182            (from, to(by))
183        } else {
184            (to(by), from)
185        };
186        let length = integrate(
187            |t| Ok(vec![speed(&longer, t)?]),
188            &[a, b],
189            1,
190            &Quadrature::default(),
191        )?
192        .value[0];
193        let off = distance.sub(length);
194        last = Some((by, length));
195        if off.could_be_equal(Scalar::ZERO) {
196            return Ok(by);
197        }
198        let step = off.div(speed(&longer, to(by))?)?;
199        let next = by.add(step).sharpen();
200        if !next.definitely_greater(Scalar::ZERO) {
201            break;
202        }
203        // A step that repeats its seed is at the iteration's fixed point,
204        // cycling in the last bit: done.
205        if next.could_be_equal(by) {
206            return Ok(by);
207        }
208        by = next;
209    }
210    Err(GeopError::new(format!(
211        "the surface does not run on {distance:?} past its edge: carried on, the last (parameter, length) was {last:?}"
212    )))
213}
214
215/// Extends the face standing on its own that the edge `edge` bounds past
216/// it by `distance` (see the module docs).
217pub fn extend(part: &mut NamedModel, edge: &str, distance: Scalar) -> GeopResult<FaceId> {
218    let ctx = with_context!("extend({edge}, {distance:?})");
219    if !distance.definitely_greater(Scalar::ZERO) {
220        return Err(GeopError::new(format!(
221            "an extension needs a distance greater than zero, not {distance:?}"
222        )))
223        .with_context(ctx);
224    }
225    let model = part.topology();
226    let edge_id = part.edge_id(edge).with_context(ctx)?;
227    let coedges = model.coedges_of_edge(edge_id);
228    let [coedge] = coedges.as_slice() else {
229        return Err(GeopError::new(format!(
230            "edge {edge} bounds {} faces: extending is past the free edge of a face standing on its own",
231            coedges.len()
232        )))
233        .with_context(ctx);
234    };
235    let face_id = model.get_coedge(*coedge)?.face;
236    let face_name = part.name(face_id).map(str::to_string)?;
237    let sheet = (part.sheet_of_face(face_id, "a sheet's face is extended")).with_context(ctx)?;
238    if model.body_faces(Body::Sheet(sheet))?.len() != 1 {
239        return Err(GeopError::new(format!(
240            "face {face_name} shares edges with other faces of its sheet: extending a face standing alone is supported"
241        )))
242        .with_context(ctx);
243    }
244    let face = model.get_face(face_id)?;
245    let not_whole = || {
246        GeopError::new(format!(
247            "face {face_name} is not bounded by its surface's four sides: extending such a face is not supported"
248        ))
249    };
250    let BoundaryType::Loop(anchor) = face.outer else {
251        return Err(not_whole()).with_context(ctx);
252    };
253    if !face.holes.is_empty() {
254        return Err(not_whole()).with_context(ctx);
255    }
256    // Per side of the patch, its coedge's edge, and per corner — where a
257    // side starts — its vertex.
258    let mut edges = [None; 4];
259    let mut corners = [None; 4];
260    for c in model.iterate_loop_coedges(anchor) {
261        let coedge = model.get_coedge(c)?;
262        let CoedgeGeometry::Edge(e) = coedge.geometry else {
263            return Err(not_whole()).with_context(ctx);
264        };
265        let Some(k) = side_of(&face.surface, &coedge.pcurve)? else {
266            let (t0, t1) = coedge.pcurve.domain();
267            return Err(not_whole().with_context(format!(
268                "edge {}: its pcurve of degree {} runs from {:?} to {:?}, on none of the sides of the domain {:?} x {:?}",
269                part.name(e).map(str::to_string)?,
270                coedge.pcurve.degree,
271                coedge.pcurve.evaluate(t0)?,
272                coedge.pcurve.evaluate(t1)?,
273                face.surface.domain_u(),
274                face.surface.domain_v()
275            )))
276            .with_context(ctx);
277        };
278        if edges[k].is_some() {
279            return Err(not_whole()).with_context(ctx);
280        }
281        edges[k] = Some(e);
282        corners[k] = Some(model.coedge_start_vertex_id(c)?);
283    }
284    let (Some(edges), Some(corners)) = (
285        edges.into_iter().collect::<Option<Vec<_>>>(),
286        corners.into_iter().collect::<Option<Vec<_>>>(),
287    ) else {
288        return Err(not_whole()).with_context(ctx);
289    };
290    let side = edges
291        .iter()
292        .position(|&e| e == edge_id)
293        .expect("the edge bounds the face");
294
295    // How far in the parameter across the edge: the distance over how fast
296    // the surface runs across it at its middle — sharpened, a free choice
297    // of how much to add.
298    let surface = &face.surface;
299    let ((u0, u1), (v0, v1)) = (surface.domain_u(), surface.domain_v());
300    let half = |a: Scalar, b: Scalar| a.add(b).div(Scalar::TWO);
301    let (along_u, at_end, (u, v)) = match side {
302        0 => (false, false, (half(u0, u1)?, v0)),
303        1 => (true, true, (u1, half(v0, v1)?)),
304        2 => (false, true, (half(u0, u1)?, v1)),
305        _ => (true, false, (u0, half(v0, v1)?)),
306    };
307    let carried = |by: Scalar| {
308        surface.extended(along_u, at_end, by).map_err(|e| {
309            e.with_context(format!(
310                "carrying the surface of face {face_name} on past edge {edge}"
311            ))
312        })
313    };
314    let by = parameter_for_length(surface, along_u, at_end, (u, v), distance, carried)
315        .with_context(ctx)?;
316    let longer = carried(by)?;
317
318    // The face anew: the whole longer patch, its sides its border curves,
319    // each keeping the name of the edge it continues.
320    let ((u0, u1), (v0, v1)) = (longer.domain_u(), longer.domain_v());
321    let border = |u_fixed: bool, first: bool| {
322        boundary_curve(&longer, u_fixed, first)
323            .ok_or_else(|| GeopError::new("the extended surface is not clamped at its border"))
324    };
325    // Sides 0..4, as the loop runs: from corner `k` to corner `k + 1`. The
326    // curves along `u` run with sides 0 and 2's `u`, along `v` with 1 and
327    // 3's `v`: 0 and 1 forward, 2 and 3 backward.
328    let curves = [
329        border(false, true)?,
330        border(true, false)?,
331        border(false, false)?,
332        border(true, true)?,
333    ];
334    let corner_points = [
335        longer.evaluate(u0, v0)?,
336        longer.evaluate(u1, v0)?,
337        longer.evaluate(u1, v1)?,
338        longer.evaluate(u0, v1)?,
339    ];
340    let uv = |u: Scalar, v: Scalar| Vector2::from_array([u, v]);
341    let corner_uv = [uv(u0, v0), uv(u1, v0), uv(u1, v1), uv(u0, v1)];
342    let spec = BodySpec {
343        vertices: corner_points.to_vec(),
344        edges: (0..4)
345            .map(|k| {
346                let forward = k < 2;
347                EdgeSpec {
348                    curve: curves[k].clone(),
349                    start: if forward { k } else { (k + 1) % 4 },
350                    end: if forward { (k + 1) % 4 } else { k },
351                }
352            })
353            .collect(),
354        faces: vec![FaceSpec {
355            surface: longer.clone(),
356            outer: (0..4)
357                .map(|k| {
358                    Ok(CoedgeSpec {
359                        on: CoedgeOn::Edge(k, Sense::from_forward(k < 2)),
360                        pcurve: NurbCurve2D::line(corner_uv[k], corner_uv[(k + 1) % 4])?,
361                    })
362                })
363                .collect::<GeopResult<_>>()?,
364            holes: Vec::new(),
365        }],
366        shells: vec![vec![0]],
367        solid: false,
368    };
369    let names = BodyNames {
370        vertices: corners
371            .iter()
372            .map(|&v| part.name(v).map(str::to_string))
373            .collect::<GeopResult<_>>()?,
374        edges: edges
375            .iter()
376            .map(|&e| part.name(e).map(str::to_string))
377            .collect::<GeopResult<_>>()?,
378        faces: vec![face_name],
379        solid: None,
380    };
381    part.delete_bodies(&[Body::Sheet(sheet)])
382        .with_context(ctx)?;
383    let built = part.build_body(spec, names).with_context(ctx)?;
384    Ok(built.faces[0])
385}
386
387#[cfg(test)]
388mod tests;
geop_ops_surface/extend.rs

geop_ops_surface/
extend.rs
