From 99d7b8c6b2f6c9925ce32807bffb256dd8e0d6a0 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Mon, 31 Mar 2025 07:47:56 +0200 Subject: [PATCH 1/8] Add extend and speed up distance calc --- cadquery/occ_impl/shapes.py | 33 ++++++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index b39efc042..eab9fea40 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3484,6 +3484,18 @@ def isolines( return [self.isoline(p, direction) for p in params] + def extend( + self, umin: bool, umax: bool, vmin: bool, vmax: bool, d: float + ) -> "Face": + """ + Extend a face. + """ + + rv = TopoDS_Face() + BRepLib.ExtendFace_s(self.wrapped, d, umin, umax, vmin, vmax, rv) + + return self.__class__(rv) + class Shell(Shape): """ @@ -5672,11 +5684,17 @@ def imprint( # collect shapes present in the history dict for k, v in history.items(): if isinstance(k, str): - history[k] = _compound_or_shape(list(images.Find(v.wrapped))) + try: + history[k] = _compound_or_shape(list(images.Find(v.wrapped))) + except Standard_NoSuchObject: + pass # store all top-level shape relations for s in shapes: - history[s] = _compound_or_shape(list(images.Find(s.wrapped))) + try: + history[s] = _compound_or_shape(list(images.Find(s.wrapped))) + except Standard_NoSuchObject: + pass return _compound_or_shape(builder.Shape()) @@ -6142,6 +6160,15 @@ def closest(s1: Shape, s2: Shape) -> Tuple[Vector, Vector]: """ Closest points between two shapes. """ - ext = BRepExtrema_DistShapeShape(s1.wrapped, s2.wrapped) + # configure + ext = BRepExtrema_DistShapeShape() + ext.SetMultiThread(True) + + # load shapes + ext.LoadS1(s1.wrapped) + ext.LoadS2(s2.wrapped) + + # perform + assert ext.Perform() return Vector(ext.PointOnShape1(1)), Vector(ext.PointOnShape2(1)) From 0e190869226fc057c528e8761c09ded0a066e8ca Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 1 Apr 2025 08:38:38 +0200 Subject: [PATCH 2/8] Add conversion to NURBS --- cadquery/occ_impl/shapes.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index eab9fea40..d73a49b1f 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3485,14 +3485,23 @@ def isolines( return [self.isoline(p, direction) for p in params] def extend( - self, umin: bool, umax: bool, vmin: bool, vmax: bool, d: float + self, d: float, umin: bool, umax: bool, vmin: bool, vmax: bool ) -> "Face": """ - Extend a face. + Extend a face. Does not work well in periodic directions. + + :param d: length of the extension. + :param umin: extend along the umin isoline. + :param umax: extend along the umax isoline. + :param vmin: extend along the vmin isoline. + :param umax: extend along the umax isoline. """ + # convert to NURBS first + tmp = self.toNURBS().wrapped + rv = TopoDS_Face() - BRepLib.ExtendFace_s(self.wrapped, d, umin, umax, vmin, vmax, rv) + BRepLib.ExtendFace_s(tmp, d, umin, umax, vmin, vmax, rv) return self.__class__(rv) From 597d63c6ef767f97f1e2991fc2eb79743fb96cab Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 1 Apr 2025 19:24:20 +0200 Subject: [PATCH 3/8] Add checks to extend --- cadquery/occ_impl/shapes.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index d73a49b1f..f1f9ee5a2 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -157,6 +157,7 @@ Geom_Plane, Geom_BSplineCurve, Geom_Curve, + Geom_BSplineSurface, ) from OCP.Geom2d import Geom2d_Line @@ -3497,11 +3498,17 @@ def extend( :param umax: extend along the umax isoline. """ - # convert to NURBS first - tmp = self.toNURBS().wrapped + # convert to NURBS if needed + tmp = self.toNURBS() if self.geomType() != "BSPLINE" else self + + # check min degree + ga = tcast(Geom_BSplineSurface, tmp._geomAdaptor()) + min_deg = min(ga.UDegree(), ga.VDegree()) + + assert min_deg >= 3, "At least degree 3 surface required" rv = TopoDS_Face() - BRepLib.ExtendFace_s(tmp, d, umin, umax, vmin, vmax, rv) + BRepLib.ExtendFace_s(tmp.wrapped, d, umin, umax, vmin, vmax, rv) return self.__class__(rv) @@ -5955,6 +5962,7 @@ def sweep( def _make_builder(): rv = BRepOffsetAPI_MakePipeShell(spine.wrapped) + if aux: rv.SetMode(_get_one_wire(aux).wrapped, True) else: From 80e2ebff28af5c7e7822c6904947f6c7aae14b22 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 1 Apr 2025 19:37:45 +0200 Subject: [PATCH 4/8] Remove checks --- cadquery/occ_impl/shapes.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index f1f9ee5a2..ea9335c5c 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3501,12 +3501,6 @@ def extend( # convert to NURBS if needed tmp = self.toNURBS() if self.geomType() != "BSPLINE" else self - # check min degree - ga = tcast(Geom_BSplineSurface, tmp._geomAdaptor()) - min_deg = min(ga.UDegree(), ga.VDegree()) - - assert min_deg >= 3, "At least degree 3 surface required" - rv = TopoDS_Face() BRepLib.ExtendFace_s(tmp.wrapped, d, umin, umax, vmin, vmax, rv) From dc9016838aa9052a7eb26b067b9801caf6b3ace4 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Tue, 1 Apr 2025 19:45:17 +0200 Subject: [PATCH 5/8] Add test and defaults --- cadquery/occ_impl/shapes.py | 7 ++++++- tests/test_shapes.py | 9 +++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index ea9335c5c..dac0c9e72 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -3486,7 +3486,12 @@ def isolines( return [self.isoline(p, direction) for p in params] def extend( - self, d: float, umin: bool, umax: bool, vmin: bool, vmax: bool + self, + d: float, + umin: bool = True, + umax: bool = True, + vmin: bool = True, + vmax: bool = True, ) -> "Face": """ Extend a face. Does not work well in periodic directions. diff --git a/tests/test_shapes.py b/tests/test_shapes.py index 5664afe78..5af2e3a26 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -13,6 +13,7 @@ cylinder, ellipse, spline, + sweep, ) from pytest import approx, raises @@ -255,3 +256,11 @@ def test_isolines(): assert isos_u[0].Length() == approx(2) assert isos_v[0].Length() == approx(pi) + + +def test_extend(): + + f = sweep(spline((0, 0), (0, 1), (2, 0)), spline((0, 0, 0), (0, 1, 1), (0, 1, 5))) + f_ext = f.extend(1) + + assert f_ext.Area() > f.Area() From 0dfc4076b713f298afef82f91d09a4d2b8a19452 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Thu, 3 Apr 2025 08:08:08 +0200 Subject: [PATCH 6/8] Better coverage --- tests/test_free_functions.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/test_free_functions.py b/tests/test_free_functions.py index d2c89cd96..117c9dd40 100644 --- a/tests/test_free_functions.py +++ b/tests/test_free_functions.py @@ -430,6 +430,7 @@ def test_imprint(): assert len(res_glue_full.Faces()) == len(compound(b1, b2).Faces()) - 1 + # imprint with history history = dict(b1=b1, b3=b3) res_glue_partial = imprint(b1, b3, glue="partial", history=history) @@ -439,6 +440,13 @@ def test_imprint(): assert len(b1_imp.Faces()) == len(b1.Faces()) + 1 assert len(res_glue_partial.Faces()) == len(b1_imp.Faces() + b3_imp.Faces()) - 1 + # imprint with faulty history + history = dict(b2=b2) + # this does not raise! + res_glue_partial = imprint(b1, b3, glue="partial", history=history) + + assert b2 not in history + def test_setThreads(): From b606d0192262415b62ea7ebff34f95f8aa79f551 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Thu, 3 Apr 2025 19:16:24 +0200 Subject: [PATCH 7/8] Ignore coverege of a failed imprint --- cadquery/occ_impl/shapes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index dac0c9e72..e72b92c99 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -5708,7 +5708,7 @@ def imprint( for s in shapes: try: history[s] = _compound_or_shape(list(images.Find(s.wrapped))) - except Standard_NoSuchObject: + except Standard_NoSuchObject: # pragma: no cover pass return _compound_or_shape(builder.Shape()) From 37f5d9359a3be3be08bbe88a46b2f67f70480987 Mon Sep 17 00:00:00 2001 From: adam-urbanczyk <13981538+adam-urbanczyk@users.noreply.github.com> Date: Fri, 4 Apr 2025 08:55:56 +0200 Subject: [PATCH 8/8] Initial commit of bsplines --- cadquery/occ_impl/bsplines.py | 277 ++++++++++++++++++++++++++++++++++ 1 file changed, 277 insertions(+) create mode 100644 cadquery/occ_impl/bsplines.py diff --git a/cadquery/occ_impl/bsplines.py b/cadquery/occ_impl/bsplines.py new file mode 100644 index 000000000..b7ff9dc11 --- /dev/null +++ b/cadquery/occ_impl/bsplines.py @@ -0,0 +1,277 @@ +from typing import Sequence + +from numpy import ( + atleast_1d, + concat, + zeros, + searchsorted, + zeros_like, +) + +from numba import njit as _njit, prange, f8 + +njit = _njit(cache=True, fastmath=False, parallel=True) +njiti = _njit(cache=True, inline="always", fastmath=False, parallel=False) + + +from cadquery.occ_impl.shapes import Edge +from cadquery.occ_impl.geom import Vector + +from OCP.Geom import Geom_BSplineCurve +from OCP.gp import gp_Pnt +from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge +from OCP.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt +from OCP.TColStd import ( + TColStd_Array1OfReal, + TColStd_Array1OfInteger, + TColStd_Array2OfPnt, +) + + +def periodicDM(knots: f8[:], xs: f8[:]) -> f8[:, :]: + """ + Periodic cubic B-spline design matrix. + """ + + order = 3 + + knots_ext = concat( + (knots[-(order + 1) : -1] - knots[-1], knots, knots[-1] + knots[1 : order + 1]) + ) + + n = len(xs) + m = len(knots) - 1 + + rv = zeros((n, m)) + + for j in prange(m + 3): + rv[:, j % m] += cubicBS(xs, j + 2, knots_ext) + + return rv, knots_ext + + +def DM(knots: f8[:], xs: f8[:]) -> f8[:, :]: + """ + Non-periodic cubic B-spline design matrix. + """ + + order = 3 + + knots_ext = concat((knots[:1].repeat(order), knots, knots[-1:].repeat(order))) + + n = len(xs) + m = len(knots) + order - 1 + + rv = zeros((n, m)) + if xs[0] == knots[0]: + rv[0, 0] = 1.0 + + for j in prange(m): + rv[:, j] += cubicBS(xs, j + order - 1, knots_ext) + + return rv, knots_ext + + +@njiti +def cubicBS(xs: f8[:], i: int, knots: f8[:]) -> f8[:]: + """ + Cubic B-spline basis function. + """ + + u = knots + ts = atleast_1d(xs) + res = zeros_like(xs) + + d1 = u[i + 1] - u[i - 2] + d2 = u[i] - u[i - 2] + d3 = u[i - 1] - u[i - 2] + d4 = u[i] - u[i - 1] + d5 = u[i + 1] - u[i - 1] + d6 = u[i + 2] - u[i - 1] + d7 = u[i + 1] - u[i] + d8 = u[i + 2] - u[i] + d9 = u[i + 2] - u[i + 1] + + j1, j2 = searchsorted(ts, (u[i - 2], u[i + 2])) + + for j in range(j1, j2 + 1): + t = ts[j] + rv = 0.0 + + if t <= u[i - 1] and t > u[i - 2]: + if d1 == 0 or d2 == 0 or d3 == 0: + rv = 0 + else: + rv = (t - u[i - 2]) ** 3 + rv /= d1 * d2 * d3 + + elif t <= u[i] and t > u[i - 1]: + tmp1 = 0.0 + + if d2 != 0 and d4 != 0: + tmp1 = (t - u[i - 2]) * (u[i] - t) + tmp1 /= d2 * d4 + + tmp2 = 0.0 + + if d5 != 0 and d4 != 0: + tmp2 = (u[i + 1] - t) * (t - u[i - 1]) + tmp2 /= d5 * d4 + + if d1 != 0: + rv += (t - u[i - 2]) / d1 + rv *= tmp1 + tmp2 + + tmp3 = 0.0 + + if d6 != 0 and d5 != 0 and d4 != 0: + tmp3 = (u[i + 2] - t) * (t - u[i - 1]) ** 2 + tmp3 /= d6 * d5 * d4 + + rv += tmp3 + + elif t <= u[i + 1] and t > u[i]: + + tmp1 = 0.0 + + if d1 != 0 and d7 != 0 and d5 != 0: + tmp1 = (t - u[i - 2]) * (u[i + 1] - t) ** 2 + tmp1 /= d1 * d7 * d5 + + rv += tmp1 + + tmp2 = 0.0 + + if d5 != 0 and d7 != 0: + tmp2 = (t - u[i - 1]) * (u[i + 1] - t) + tmp2 /= d5 * d7 + + tmp3 = 0 + + if d8 != 0 and d7 != 0: + tmp3 = (u[i + 2] - t) * (t - u[i]) + tmp3 /= d8 * d7 + + if d6 != 0: + rv += (tmp2 + tmp3) * (u[i + 2] - t) / d6 + + elif t <= u[i + 2] and t > u[i + 1]: + + if d9 != 0 and d8 != 0 and d6 != 0: + rv += (u[i + 2] - t) ** 3 + rv /= d9 * d8 * d6 + + res[j] = rv + + return res + + +def _to_pts_arr(pts: Sequence[gp_Pnt]) -> TColgp_Array1OfPnt: + """ + Helper to construct a TColgp_Array1OfPnt. + """ + + rv = TColgp_Array1OfPnt(1, len(pts)) + + for i, p in enumerate(pts): + rv.SetValue(i + 1, gp_Pnt(*p)) + + return rv + + +def _to_pts_arr2(pts: Sequence[Sequence[Vector]]) -> TColStd_Array2OfPnt: + """ + Helper to construct a TColgp_Array2OfReal. + """ + + rv = TColgp_Array2OfPnt(1, len(pts), 1, len(pts[0])) + + for i, ps_v in enumerate(pts): + for j, pt in enumerate(ps_v): + rv.SetValue(i + 1, j + 1, gp_Pnt(*pt)) + + return rv + + +def _to_real_arr(vals: Sequence[float]) -> TColStd_Array1OfReal: + """ + Helper to construct a TColStd_Array1OfReal. + """ + + rv = TColStd_Array1OfReal(1, len(vals)) + + for i, v in enumerate(vals): + rv.SetValue(i + 1, v) + + return rv + + +def _nominal_mults_periodic(knots: Sequence[float]) -> TColStd_Array1OfInteger: + """ + Generate a nominal multipicity vector for a given periodic knot vector. + """ + + n_knots = len(knots) + + # init multiplicites + mults = TColStd_Array1OfInteger(1, n_knots) + mults.Init(1) + + return mults + + +def _nominal_mults(knots: Sequence[float], order: int) -> TColStd_Array1OfInteger: + """ + Generate a nominal multipicity vector for a given clamped knot vector. + """ + + n_knots = len(knots) + + # init multiplicites + mults = TColStd_Array1OfInteger(1, n_knots) + mults.Init(1) + + mults.SetValue(1, order + 1) + mults.SetValue(n_knots, order + 1) + + return mults + + +def toPeriodicCurve( + pts: Sequence[Sequence[float]], knots: Sequence[float], order: int = 3 +) -> Geom_BSplineCurve: + """ + Construct a Geom periodic curve. + """ + + return Geom_BSplineCurve( + _to_pts_arr(pts), + _to_real_arr(knots), + _nominal_mults_periodic(knots), + order, + True, + ) + + +def toCurve( + pts: Sequence[gp_Pnt], knots: Sequence[float], order: int = 3 +) -> Geom_BSplineCurve: + """ + Construct a Geom non-periodic curve. + """ + + return Geom_BSplineCurve( + _to_pts_arr(pts), + _to_real_arr(knots), + _nominal_mults(knots, order), + order, + False, + ) + + +def toEdge(geom: Geom_BSplineCurve) -> Edge: + """ + Construct an edge from Geom curve + """ + + return Edge(BRepBuilderAPI_MakeEdge(geom).Shape())