diff --git a/cadquery/func.py b/cadquery/func.py index d1fc8450e..ef65d674d 100644 --- a/cadquery/func.py +++ b/cadquery/func.py @@ -44,6 +44,7 @@ extrude, revolve, offset, + offset2D, sweep, loft, check, diff --git a/cadquery/occ_impl/nurbs.py b/cadquery/occ_impl/nurbs.py index 8386062a5..4a958d35f 100644 --- a/cadquery/occ_impl/nurbs.py +++ b/cadquery/occ_impl/nurbs.py @@ -11,7 +11,7 @@ from numpy.typing import NDArray from numpy import linspace, ndarray -from casadi import ldl, ldl_solve +from casadi import ldl, ldl_solve, Linsol, DM, solve from OCP.Geom import Geom_BSplineCurve, Geom_BSplineSurface from OCP.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt @@ -1454,7 +1454,9 @@ def uIsoMatrix(surf: Surface, u: float) -> COO: nu = surf.pts.shape[0] nv = surf.pts.shape[1] - block = designMatrix(np.atleast_1d(u), surf.uorder, surf.uknots, surf.uperiodic) + block = designMatrix( + np.atleast_1d(np.array(u)), surf.uorder, surf.uknots, surf.uperiodic + ) shape = (nv, nu * nv) @@ -1481,7 +1483,9 @@ def vIsoMatrix(surf: Surface, v: float) -> COO: nu = surf.pts.shape[0] nv = surf.pts.shape[1] - block = designMatrix(np.atleast_1d(v), surf.vorder, surf.vknots, surf.vperiodic) + block = designMatrix( + np.atleast_1d(np.array(v)), surf.vorder, surf.vknots, surf.vperiodic + ) shape = (nu, nu * nv) @@ -1853,8 +1857,9 @@ def approximate2D( CtC = C.T @ C # solve normal equations - D, L, P = ldl(CtC, False) - pts = ldl_solve(C.T @ data, D, L, P).toarray() + CtC = DM(CtC) + ldl = Linsol("mumps", "mumps", CtC.sparsity()) + pts = ldl.solve(CtC, C.T @ data).toarray() # construct the result rv = Surface( @@ -1927,8 +1932,10 @@ def constrainedApproximate2D( LHS = sp.bmat(((CtC_xyz, As.coo().T,), (As.coo(), None))) RHS = np.concatenate((C_xyz.T @ data.flatten(order="F"), bs)) - D, L, P = ldl(LHS, False) - sol = ldl_solve(RHS, D, L, P).toarray() + LHS_DM = DM(LHS) + ldl = Linsol("mumps", "mumps", LHS_DM.sparsity()) + sol = ldl.solve(LHS_DM, RHS).toarray() + pts = sol[: CtC_xyz.shape[0]] # construct the result @@ -2006,8 +2013,9 @@ def constrainedApproximate2D( LHS = CtC RHS = C.T @ data[:, i] - D, L, P = ldl(LHS, False) - sol = ldl_solve(RHS, D, L, P).toarray() + LHS_DM = DM(LHS) + ldl = Linsol("mumps", "mumps", LHS_DM.sparsity()) + sol = ldl.solve(LHS_DM, RHS).toarray() pts.append(sol[: CtC.shape[0]]) # construct the result diff --git a/cadquery/occ_impl/shapes.py b/cadquery/occ_impl/shapes.py index 6055b08c3..38b641372 100644 --- a/cadquery/occ_impl/shapes.py +++ b/cadquery/occ_impl/shapes.py @@ -21,10 +21,13 @@ from io import BytesIO +from warnings import warn from vtkmodules.vtkCommonDataModel import vtkPolyData from vtkmodules.vtkFiltersCore import vtkTriangleFilter, vtkPolyDataNormals +from OCP.ShapeBuild import ShapeBuild_ReShape + from .geom import Vector, VectorLike, BoundBox, Plane, Location, Matrix from .shape_protocols import geom_LUT_FACE, geom_LUT_EDGE, Shapes, Geoms @@ -1734,7 +1737,6 @@ def _repr_javascript_(self): def __iter__(self) -> Iterator["Shape"]: """ Iterate over subshapes. - """ it = TopoDS_Iterator(self.wrapped) @@ -1743,31 +1745,31 @@ def __iter__(self) -> Iterator["Shape"]: yield Shape.cast(it.Value()) it.Next() - def ancestors(self, shape: "Shape", kind: Shapes) -> "Compound": + def ancestors(self, ctx: "Shape", kind: Shapes) -> "Compound": """ - Iterate over ancestors, i.e. shapes of same kind within shape that contain self. - + Iterate over ancestors, i.e. shapes of type kind within ctx shape that contain self. """ shape_map = TopTools_IndexedDataMapOfShapeListOfShape() TopExp.MapShapesAndAncestors_s( - shape.wrapped, shapetype(self.wrapped), inverse_shape_LUT[kind], shape_map + ctx.wrapped, shapetype(self.wrapped), inverse_shape_LUT[kind], shape_map ) return Compound.makeCompound( Shape.cast(s) for s in shape_map.FindFromKey(self.wrapped) ) - def siblings(self, shape: "Shape", kind: Shapes, level: int = 1) -> "Compound": + def siblings( + self, ctx: "Shape", kind: Shapes, level: int | Iterable[int] = 1 + ) -> "Compound": """ - Iterate over siblings, i.e. shapes within shape that share subshapes of kind with self. - + Iterate over siblings, i.e. shapes within ctx shape that share subshapes of type kind with self. """ shape_map = TopTools_IndexedDataMapOfShapeListOfShape() TopExp.MapShapesAndAncestors_s( - shape.wrapped, inverse_shape_LUT[kind], shapetype(self.wrapped), shape_map, + ctx.wrapped, inverse_shape_LUT[kind], shapetype(self.wrapped), shape_map, ) exclude = TopTools_MapOfShape() @@ -1789,7 +1791,17 @@ def _siblings(shapes, level): return rv if level == 1 else _siblings(rv, level - 1) - return Compound.makeCompound(_siblings([self], level)) + # simple query + if isinstance(level, int): + rv = _siblings([self], level) + else: + # collect all the requested levels + rv = set() + for lvl in level: + exclude.Clear() + rv |= _siblings([self], lvl) + + return Compound.makeCompound(rv) def __add__(self, other: "Shape") -> "Shape": """ @@ -1933,6 +1945,20 @@ def __getitem__(self, item: int | slice) -> "Shape": else: return list(self)[item] + def size(self) -> int: + """ + Simple size implementation. + """ + + return len(list(self)) + + def reverse(self) -> "Shape": + """ + Return a copy of self with reversed orientation. + """ + + return self.cast(self.wrapped.Reversed()) + class ShapeProtocol(Protocol): @property @@ -4873,9 +4899,9 @@ def intersect( return tcast(Compound, self._bool_op(self, toIntersect, intersect_op)) - def ancestors(self, shape: "Shape", kind: Shapes) -> "Compound": + def ancestors(self, ctx: "Shape", kind: Shapes) -> "Compound": """ - Iterate over ancestors, i.e. shapes of same kind within shape that contain elements of self. + Iterate over ancestors, i.e. shapes of same kind within ctx shape that contain elements of self. """ @@ -4884,14 +4910,16 @@ def ancestors(self, shape: "Shape", kind: Shapes) -> "Compound": for t in shapetypes: TopExp.MapShapesAndAncestors_s( - shape.wrapped, t, inverse_shape_LUT[kind], shape_map + ctx.wrapped, t, inverse_shape_LUT[kind], shape_map ) return Compound.makeCompound( Shape.cast(a) for s in self for a in shape_map.FindFromKey(s.wrapped) ) - def siblings(self, shape: "Shape", kind: Shapes, level: int = 1) -> "Compound": + def siblings( + self, ctx: "Shape", kind: Shapes, level: int | Iterable[int] = 1 + ) -> "Compound": """ Iterate over siblings, i.e. shapes within shape that share subshapes of kind with the elements of self. @@ -4902,7 +4930,7 @@ def siblings(self, shape: "Shape", kind: Shapes, level: int = 1) -> "Compound": for t in shapetypes: TopExp.MapShapesAndAncestors_s( - shape.wrapped, inverse_shape_LUT[kind], t, shape_map, + ctx.wrapped, inverse_shape_LUT[kind], t, shape_map, ) exclude = TopTools_MapOfShape() @@ -4924,7 +4952,17 @@ def _siblings(shapes, level): return rv if level == 1 else _siblings(rv, level - 1) - return Compound.makeCompound(_siblings(self, level)) + # simple query + if isinstance(level, int): + rv = _siblings(self, level) + else: + # collect all the requested levels + rv = set() + for lvl in level: + exclude.Clear() + rv |= _siblings(self, lvl) + + return Compound.makeCompound(rv) def sortWiresByBuildOrder(wireList: List[Wire]) -> List[List[Wire]]: @@ -5276,6 +5314,44 @@ def _shape(s: TopoDS_Shape, _: type[T]) -> T: return tcast(T, Shape.cast(s)) +def _shape_to_faces_shells( + s: TopoDS_Shape, +) -> tuple[list[TopoDS_Face], list[TopoDS_Shell]]: + """ + Split a TopoDS_Shape into Faces and Shells. Raise if another type was found. + """ + + rv_faces = [] + rv_shells = [] + + t = s.ShapeType() + + if t == TopAbs_ShapeEnum.TopAbs_COMPOUND: + it = TopoDS_Iterator(s) + + while it.More(): + el = it.Value() + + t = el.ShapeType() + + if t == TopAbs_ShapeEnum.TopAbs_FACE: + rv_faces.append(TopoDS.Face(el)) + elif t == TopAbs_ShapeEnum.TopAbs_SHELL: + rv_shells.append(TopoDS.Shell(el)) + else: + raise ValueError(f"Found unexpected type {t}") + + it.Next() + elif t == TopAbs_ShapeEnum.TopAbs_SHELL: + rv_shells.append(TopoDS.Shell(s)) + elif t == TopAbs_ShapeEnum.TopAbs_FACE: + rv_faces.append(TopoDS.Face(s)) + else: + raise ValueError(f"Found unexpected type {t}") + + return rv_faces, rv_shells + + def _pts_to_harray(pts: Sequence[VectorLike]) -> TColgp_HArray1OfPnt: """ Convert a sequence of Vector to a TColgp harray (OCCT specific). @@ -5595,7 +5671,7 @@ def shell( manifold: bool = True, ctx: Optional[Sequence[Shape] | Shape] = None, history: Optional[ShapeHistory] = None, -) -> Shell: +) -> Shape: """ Build shell from faces. If ctx is specified, local sewing is performed. """ @@ -5618,20 +5694,24 @@ def shell( sewed = builder.SewedShape() _process_sewing_history(builder, faces, history) - rv: Union[TopoDS_Shape, TopoDS_Shell] + rv = [] - # for one face sewing will not produce a shell - if sewed.ShapeType() == TopAbs_ShapeEnum.TopAbs_FACE: - rv = TopoDS_Shell() + # split the results + rv_faces, rv_shells = _shape_to_faces_shells(sewed) + + rv.extend(rv_shells) + + # for one closed face sewing will not produce a shell, so a special treatment is needed + for f_topo in rv_faces: + tmp = TopoDS_Shell() builder_topo = TopoDS_Builder() - builder_topo.MakeShell(rv) - builder_topo.Add(rv, sewed) + builder_topo.MakeShell(tmp) + builder_topo.Add(tmp, f_topo) - else: - rv = sewed + rv.append(tmp) - return _shape(rv, Shell) + return _compound_or_shape(rv) @multidispatch @@ -5641,7 +5721,7 @@ def shell( manifold: bool = True, ctx: Optional[Sequence[Shape] | Shape] = None, history: Optional[ShapeHistory] = None, -) -> Shell: +) -> Shape: """ Build shell from a sequence of faces. If ctx is specified, local sewing is performed. """ @@ -5693,9 +5773,17 @@ def solid( builder.Add(sh.wrapped) # fix orientations + ctx = ShapeBuild_ReShape() + sf = ShapeFix_Solid(builder.Solid()) + sf.SetContext(ctx) sf.Perform() + # update history if applicable + if history is not None: + for k, v in history.items(): + history[k] = Shape.cast(ctx.Apply(v.wrapped)) + return _shape(sf.Solid(), Solid) @@ -6459,6 +6547,48 @@ def _offset(t): return rv +def offset2D( + s: Shape, + t: float, + ctx: Shape | None = None, + kind: Literal["arc", "intersection", "tangent"] = "arc", + closed: bool = True, +) -> Shape: + """ + 2D offset edges or wires. ctx face might be needed for ambiguous wires/edges. + Only works with planar geometries. + """ + + kind_dict = { + "arc": GeomAbs_JoinType.GeomAbs_Arc, + "intersection": GeomAbs_JoinType.GeomAbs_Intersection, + "tangent": GeomAbs_JoinType.GeomAbs_Tangent, + } + + if ctx: + # build a dummy face based on the geometry of the ctx face. + fbldr = BRepBuilderAPI_MakeFace(_get_one(ctx, "Face")._geomAdaptor(), 1e-6) + + for el in _get_wires(s): + fbldr.Add(el.wrapped) + + fbldr.Build() + + bldr = BRepOffsetAPI_MakeOffset(fbldr.Face(), kind_dict[kind], not closed) + + else: + # simple case - input wire defines a plane + bldr = BRepOffsetAPI_MakeOffset() + bldr.Init(kind_dict[kind], not closed) + + for el in _get_wires(s): + bldr.AddWire(el.wrapped) + + bldr.Perform(t) + + return _compound_or_shape(bldr.Shape()) + + @multimethod def sweep( s: Shape, path: Shape, aux: Optional[Shape] = None, cap: bool = False @@ -6748,13 +6878,14 @@ def project( return _normalize(compound(results)) -#%% diagnotics +#%% diagnostics def check( s: Shape, results: Optional[List[Tuple[List[Shape], Any]]] = None, tol: Optional[float] = None, + verbose: bool = True, ) -> bool: """ Check if a shape is valid. @@ -6771,14 +6902,19 @@ def check( if tol: analyzer.SetFuzzyValue(tol) - # output detailed results if requested + # generate warnings, output detailed results if requested if results is not None: results.clear() - for r in analyzer.Result(): - results.append( - (_toptools_list_to_shapes(r.GetFaultyShapes1()), r.GetCheckStatus()) - ) + for r in analyzer.Result(): + res = (_toptools_list_to_shapes(r.GetFaultyShapes1()), r.GetCheckStatus()) + msg = f"\n\tCheck failed.\n\tSubshapes: {res[0]} \n\tStatus: {res[1]}" + + if verbose: + warn(msg) + + if results is not None: + results.append(res) return rv diff --git a/tests/test_free_functions.py b/tests/test_free_functions.py index afd07a9a7..d63711419 100644 --- a/tests/test_free_functions.py +++ b/tests/test_free_functions.py @@ -46,6 +46,7 @@ project, edgeOn, faceOn, + offset2D, ) from cadquery.occ_impl.shapes import ( @@ -55,6 +56,7 @@ _get_one, _get_edges, _adaptor_curve_to_edge, + _shape_to_faces_shells, ) from OCP.BOPAlgo import BOPAlgo_CheckStatus @@ -143,6 +145,16 @@ def test_adaptor_curve_to_edge(): assert isinstance(e, TopoDS_Edge) +def test__shape_to_faces_shells(): + + # bad weather tests + with raises(ValueError): + _shape_to_faces_shells(compound(vertex(1, 0, 0), vertex(0, 0, 0)).wrapped) + + with raises(ValueError): + _shape_to_faces_shells(vertex(0, 0, 0).wrapped) + + # %% constructors @@ -233,6 +245,10 @@ def test_sewing(): def test_solid(): b = box(1, 1, 1) + b_large = box(10, 10, 1) + b_small = box(0.1, 0.1, 0.1).moved(b_large) + sphere1 = sphere(0.1).moved(b_large) + sphere2 = sphere1.moved(x=2) # solid s1 = solid(b.Faces()) @@ -264,6 +280,18 @@ def test_solid(): for f in final_faces: assert f in final_faces_history + # solid with multiple periodic voids + s5 = solid(b_large.Faces(), inner=sphere1.Faces() + sphere2.Faces()) + + assert s5.isValid() + assert s5.Volume() == approx(b_large.Volume() - 2 * sphere1.Volume()) + + # solid with multiple simple voids + s6 = solid(b_large.Faces(), inner=b_small.Faces() + b_small.moved(x=1).Faces()) + + assert s6.isValid() + assert s6.Volume() == approx(b_large.Volume() - 2 * b_small.Volume()) + def test_edgeOn(): @@ -801,6 +829,34 @@ def test_offset(): assert r4.Volume() == approx(4) +def test_offset2D(): + + c = circle(1) + seg = segment((0, 0), (1, 1)) + ctx = plane() + + # simple offset + r1 = offset2D(c, 0.1) + + # offset with a context + r2 = offset2D(seg, 0.2, ctx) + + # offset that is not closed + r3 = offset2D(seg, -0.3, ctx, closed=False) + + assert r1.isValid() + assert len(clean(r1).Edges()) == 1 + assert face(r1).isValid() + + assert r2.isValid() + assert len(clean(r2).Edges()) == 4 + assert face(r2).isValid() + + assert r3.isValid() + assert len(clean(r3).Edges()) == 1 + assert r3.edge().Length() == approx(seg.Length()) + + def test_sweep(): w1 = rect(1, 1) diff --git a/tests/test_nurbs.py b/tests/test_nurbs.py index 277608370..a3255ef66 100644 --- a/tests/test_nurbs.py +++ b/tests/test_nurbs.py @@ -586,8 +586,8 @@ def test_constrainedApproximate2D(torus_surf, lam, penalty): assert not np.allclose(isov.pts[:, 2], 0) # # constraints per direction - Au = uIsoMatrix(surf, np.array(0.0)) - Av = vIsoMatrix(surf, np.array(0.0)) + Au = uIsoMatrix(surf, 0.0) + Av = vIsoMatrix(surf, 0.0) by = np.zeros(Au.shape[0]) bz = np.zeros(Av.shape[0]) diff --git a/tests/test_shapes.py b/tests/test_shapes.py index 006a871e2..4fae3f107 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -23,15 +23,22 @@ polygon, wireOn, vertex, + fuse, ) from cadquery.selectors import NearestToPointSelector -from pytest import approx, raises +from pytest import approx, raises, fixture from math import pi +@fixture +def simple_box(): + + return box(1, 1, 1) + + def test_edge_paramAt(): # paramAt for a segment @@ -404,3 +411,60 @@ def test_center(): assert c.x == approx(1) assert c.y == approx(1) assert c.z == approx(1) + + +def test_reverse(simple_box): + + simple_box = box(1, 1, 1) + + assert simple_box.Volume() > 0 + + br = simple_box.reverse() + + # reversed solid will have a negative volume + assert br.Volume() < 0 + + # normals will be pointing in opposite direction + delta = simple_box.face().normalAt() + br.face().normalAt() + assert delta.Length == approx(0) + + +def test_siblings(simple_box): + + f = simple_box.face("Z").Center()).Length == approx(0) + + siblings_cmp_12 = simple_box.faces(">Z").siblings(simple_box, "Edge", (1, 2)) + assert siblings_cmp_12.size() == 5 + + siblings_cmp_2 = simple_box.faces(">Z").siblings(simple_box, "Edge", (2,)) + assert siblings_cmp_2.size() == 1 + + siblings_cmp_edges_12 = simple_box.edges(">Z").siblings( + simple_box, "Vertex", (1, 2) + ) + assert siblings_cmp_edges_12.size() == 8 + + # more complex shape + stacked_box = fuse( + simple_box, simple_box.moved(x=1), simple_box.moved(x=2), simple_box.moved(x=3), + ) + + f = stacked_box.face("