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()) diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index b39efc042..e72b92c99 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 @@ -3484,6 +3485,32 @@ def isolines( return [self.isoline(p, direction) for p in params] + def extend( + 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. + + :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 if needed + tmp = self.toNURBS() if self.geomType() != "BSPLINE" else self + + rv = TopoDS_Face() + BRepLib.ExtendFace_s(tmp.wrapped, d, umin, umax, vmin, vmax, rv) + + return self.__class__(rv) + class Shell(Shape): """ @@ -5672,11 +5699,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: # pragma: no cover + pass return _compound_or_shape(builder.Shape()) @@ -5928,6 +5961,7 @@ def sweep( def _make_builder(): rv = BRepOffsetAPI_MakePipeShell(spine.wrapped) + if aux: rv.SetMode(_get_one_wire(aux).wrapped, True) else: @@ -6142,6 +6176,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)) 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(): 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()