Assembly support

[adam-urbanczyk]

This PR will eventually resolve #276, #20

What's left

• Exact jacobian , finite diff jacobian by hand
• Finalize solver choice -> BFGS it is
• Implement parameters handling (e.g. specific angle or distance constraint)
• Visually attractive example in docs/README

[codecov]

Codecov Report

Merging #440 into master will increase coverage by 0.25%.
The diff coverage is 93.20%.

@@            Coverage Diff             @@
##           master     #440      +/-   ##
==========================================
+ Coverage   93.61%   93.87%   +0.25%     
==========================================
  Files          25       30       +5     
  Lines        5326     6984    +1658     
  Branches      554      807     +253     
==========================================
+ Hits         4986     6556    +1570     
- Misses        215      259      +44     
- Partials      125      169      +44     

Impacted Files	Coverage Δ
cadquery/occ_impl/geom.py	88.93% <75.00%> (+0.09%)	⬆️
cadquery/assembly.py	88.43% <88.43%> (ø)
cadquery/occ_impl/solver.py	89.90% <89.90%> (ø)
cadquery/occ_impl/assembly.py	95.06% <95.06%> (ø)
cadquery/__init__.py	100.00% <100.00%> (ø)
cadquery/occ_impl/exporters/assembly.py	100.00% <100.00%> (ø)
cadquery/occ_impl/shapes.py	90.89% <100.00%> (+0.14%)	⬆️
tests/test_assembly.py	100.00% <100.00%> (ø)
tests/test_cad_objects.py	98.99% <100.00%> (+0.02%)	⬆️
tests/test_cadquery.py	98.83% <100.00%> (-0.15%)	⬇️
... and 8 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2d721d0...21f7f48. Read the comment docs.

[jmwright]

Tagging @bernhard-42 on this since jupyter-cadquery supports lightweight assemblies, and it may be useful to compare notes.

[bernhard-42]

When I look at the above examples, cq.Location(cq.Vector(x, y, z)) appears quite often. Would it be possible to support short cuts like

assy3 = cq.Assembly(b3, loc=(0, 0, 5), name="THIRD", color='orange')

and let the __init__ function create the correct Location and Color object.

[bernhard-42]

This is how I would visualize the example in jupyter-cadquery:
I see one TOP assembly and 3 sub assemblies (SECOND, THIRD and 4th). And every assembly can have several objects (object_0, object_1, ...):

From how it is built I would have seen a difference between SECOND/THIRD and 4th: SECOND and THIRD are assemblies itself, while 4th is a shape. I would have expected 4th to be on the same level as the shapes of TOP.
So it looks like

assy = cq.Assembly(b1, loc=L(1, 1, 0), name="TOP")
assy2 = cq.Assembly(b2, loc=L(0, -3, 4), name="SECOND")
assy3 = cq.Assembly(b3, loc=L(0, 0, 5), name="THIRD", color=C('orange'))

assy.add(assy2, color=C("green"))
assy.add(assy3)
assy.add(b4, loc=L(0, 0, 5), name="4th",color=C("blue1"))

assy = cq.Assembly(b1, loc=L(1, 1, 0), name="TOP")
assy2 = cq.Assembly(b2, loc=L(0, -3, 4), name="SECOND")
assy3 = cq.Assembly(b3, loc=L(0, 0, 5), name="THIRD", color=C('orange'))
assy4 = cq.Assembly(b4, loc=L(0, 0, 5), name="4th",color=C("blue1"))

assy.add(assy2, color=C("green"))
assy.add(assy3)
assy.add(assy4)

lead to identical results, which I found a little bit confusing (btw. C and L are just shortcut functions).

[bernhard-42]

On a completely different topic, I think I have now stared for two hours into the one example above, commenting out constraints and visualising constraints and results - and I must admit, I don't understand it. Constraint based assembly is totally new to me. Is there any documentation around this?

[bernhard-42]

btw. converting from cadquery assemblies to jupyter-cadquery assemblies is straightforward (however, one needs to use the Assembly and Part classes from the occ version of juypter-cadquery - need to check how to properly integrate later)

from jupyter_cadquery.occ import Part, Assembly

def rgb(assy):
    def b(x):
        return int(255*x)
    
    if assy.color is None:
        return "#aaa"
    rgb = assy.color.wrapped.GetRGB()
    return "#%02x%02x%02x" % (b(rgb.Red()), b(rgb.Green()), b(rgb.Blue()))

def convert(assy, loc=None):
    loc = assy.loc if loc is None else loc * assy.loc
    color = rgb(assy)
    parent = [Part(shape.located(loc).wrapped, "object_%d" % (i), color=color) for i, shape in enumerate(assy.shapes)]
    children = [convert(c, loc) for c in assy.children]
    return Assembly(parent + children, assy.name)

@adam-urbanczyk Does the convert function look right for you?

[adam-urbanczyk]

Thanks @bernhard-42 , looks cool! You are right regarding the product of locations. The semantics of the Assembly class is essentially contained in this snippet:

class Assembly(object):
    loc: Location
    name: str
    color: Optional[Color]
    metadata: Dict[str, Any]

    obj: Union[Shape, Workplane, None]
    parent: Optional["Assembly"]
    children: List["Assembly"]

So every child of an Assembly is an Assembly object. It may or may not contain a cq.Workplane or cq.Shape. In the view of this THIRD and 4th are actually equivalent.

Regarding constraints, this is actually a standard feature of CAD software. I have not seen it yet in scripted CAD packages (not that there are so many). For now it is highly experimental, may break and change at any time kind of thing. So no docs at the moment.

Convert looks good to me.

[bernhard-42]

Thanks, Adam

[bernhard-42]

I did my homework about constraint based assemblies looking at onshape and the Assembly3 FreeCad plugin.
I think I get the idea, so I tried to build something and selected, inspired by cqparts, a bearing:

def ring(inner_radius, outer_radius, width):
    ring = (cq.Workplane(origin=(0, 0, -width / 2))
        .circle(outer_radius).circle(inner_radius)
        .extrude(width)
    )
    return ring

tol = 0.05
ball_diam = 5

r1, r2, r3, r4 = 4, 6, 8, 10
r5 = (r3 + r2) / 2
inner_ring = ring(r1, r2, ball_diam)
outer_ring = ring(r3, r4, ball_diam)

torus = cq.CQ(cq.Solid.makeTorus(r5, ball_diam / 2 + tol))
ball = cq.Workplane().sphere(ball_diam / 2)

inner = inner_ring.cut(torus)
outer = outer_ring.cut(torus)

For the issue I had with this, see #458

After having done all the steps below, I am not sure this is the best example for the assembly feature. However, you asked for some feedback on the API and so I thought to share some of my experiences:

• The API is clear and I like that it reuses the selector syntax. For the constraints I think I personally would prefer ("qos" meaning "query_or_shape"):

  def constrain(
      self, 
      id1: str,
      qos1: Union[str, Shape], 
      id2: str,
      qos2: Union[str, Shape], 
      kind: ConstraintKinds, 
      param: Any = None
  ) -> "Assembly"  

  allowing to mix the two methods

  constrain("name1", "faces@>Z", "name2", "faces@>Y", kind)
  constrain("name1", "faces@>Z", "name2", shape2, kind)
  constrain("name1", shape1, "name2", "faces@>Y", kind)
  constrain("name1", shape1, "name2", shape2, kind)

  and to me faces(">Y") and "faces@>Y" would be a better visual match than "name@faces@>Y"

• With the query syntax I haven't found a way to use points that are no vertices (e.g. center of a shape, middle points of an edge, ...). But that's maybe due to my lack of knowledge around selectors.

• The "Plane" constraint seems to always place the two faces onto each other with the shapes being stacked. How would I assemble two shapes with both faces in one plane sharing the center of the faces and having their normals showing into the same direction? I had to introduce a helper object to achieve this.

• My last example converged (after increasing maxiter in the code), however the result is not exactly correct. What is the best way to avoid this?

Finally, I don't know whether the below approach using helper objects in assemblies is appropriate, but at least it worked. And I am happy to learn, how to avoid them if possible.

The Bearing assembly

For the way I wanted to create the constraints, I needed a way to create a constraint for the center of a shape and to
get all results of an assembly query. So I built two helpers:

def center(assy, name):
    return cq.Vertex.makeVertex(*assy.objects[name].obj.val().Center().toTuple())

def query_all(assy, q):
    name, kind, arg = q.split("@")

    tmp = cq.Workplane()
    obj = assy.objects[name].obj

    if isinstance(obj, (cq.Workplane, cq.Shape)):
        tmp.add(obj)
        res = getattr(tmp, kind)(arg)

    return res.objects if isinstance(res.val(), cq.Shape) else None

Given this, I could create the assembly and add the constraints by additionally adding two construction objects:

• _center: a cylinder to use with a "Plane" constraint to arrange the rings concentric on one plane (I didn't find another way to assemble two shapes with their faces on the same plane showing into the same direction and sharing the center of the faces).
• _points: a polygon on which upper vertices the ball centers would be placed.

def balls(i):
    return "ball_%d" % i

number_balls = 5

# components
assy = cq.Assembly(outer, loc=L(0, 0, ball_diam/2), name="outer", color=C("orange"))
assy.add(inner, loc=L(20, 0, 0), name="inner", color=C("orange"))
for i in range(number_balls):
    assy.add(ball, loc=L(6*i, 20, 0), name=balls(i), color=C("black"))

# helpers
assy.add(cq.Workplane().circle(1).extrude(1), loc=L(0,-20,0), name="_center")
assy.add(cq.Workplane().polygon(number_balls, 2*r5).extrude(ball_diam/2), loc=L(20,-20,0), name="_points")

# contraints
points = query_all(assy, "_points@vertices@>Z")
cs = [
    ("outer@faces@<Z", "_center@faces@>Z", "Plane"),
    ("inner@faces@<Z", "_center@faces@>Z", "Plane"),
    ("_points@faces@<Z", "_center@faces@>Z", "Plane"),
] + [("_points", points[i],  
      balls(i), center(assy, balls(i)),  
      "Point") 
     for i in range(number_balls)]

for c in cs:
    assy.constrain(*c)

This is the assembly before solving:

And this is the result after solving:

However, when I added a 6th ball, the optimzation algorithm did not finish after after 1000 iterations. So I increased the max iterations to 10000. It finally converged

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
   60   5274   5814      2     0     0   7.478D-05   2.777D-05
  F =  2.77707271466398639E-005

CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

but the result was not exact:

[bernhard-42]

Even for 5 balls, the result is not exact:

The assembled center coordinates (after solving):

loc = assy.loc
for i in range(number_balls):
    o = assy.objects[balls(i)]
    v = o.shapes[0].located(loc * o.loc).Center()
    print("%7.4f, %7.4f, %7.4f" % (v.x, v.y, v.z))

#  6.9230, -0.9758,  2.5239
#  3.0724,  6.2966,  2.4826
# -5.0340,  4.8817,  2.4654
# -6.1935, -3.2651,  2.4959
#  1.1963, -6.8853,  2.5321

The coordinates of the polygon vertices (also after solving)

["%7.4f, %7.4f, %7.4f" % (v.X, v.Y, v.Z) for v in query_all(assy, "_points@vertices@>Z")]
# [' 7.0000,  0.0000,  2.5000',
#  ' 2.1631,  6.6574,  2.5000',
#  '-5.6631,  4.1145,  2.5000',
#  '-5.6631, -4.1145,  2.5000',
#  ' 2.1631, -6.6574,  2.5000']

I had expected that the "Point" constraint ensures that the points are colocated.
Did I dow something wrong with my constraint?

[adam-urbanczyk]

Thanks for the suggestions regarding the API. I think we can add those in future iterations of this functionality.

I changed the cost function slightly @bernhard-42 , seems to give better result now. Pleas pull and try. In general I use a numerical solver here so nothing is "exact".

BTW: please share complete code examples so that I can reproduce next time, it took me a while to glue the snippets together.

[bernhard-42]

@adam-urbanczyk The results are much better now. The issue before the change was that the solver did not align the normals well enough. To compare I created a little metrics function. For the 6 ball bearing I get

Constraint                                 Normal-Dist    Point-Dist  
--------------------------------------------------------------------
Plane outer@faces@<Z - _center@faces@>Z      0.0001571     0.0000001
Plane inner@faces@<Z - _center@faces@>Z      0.0001115     0.0000000
Plane _points@faces@<Z - _center@faces@>Z    0.0000002     0.0000001
Point _points<Vertex> - ball_0<Vertex>                     0.0000002
Point _points<Vertex> - ball_1<Vertex>                     0.0000001
Point _points<Vertex> - ball_2<Vertex>                     0.0000001
Point _points<Vertex> - ball_3<Vertex>                     0.0000002
Point _points<Vertex> - ball_4<Vertex>                     0.0000002
Point _points<Vertex> - ball_5<Vertex>                     0.0000001

with Point-Dist being norm(center1-center2) and Normal-Dist being norm(normal1 - (-normal2)) after solving.

Before Normal-Dist was up to 0.017, so for my tests the result is more than factor 100 better now - and the algorithm is much faster for this example.

All code including the metrics function can be found in https://github.com/bernhard-42/jupyter-cadquery/blob/master/examples/cadquery-assemblies.ipynb

[adam-urbanczyk]

Thanks @bernhard-42 . FYI I implemented parameter handling for Axis and Point - you can now specify the distance (default is 0) or angle (default is pi). Note that it can lead to non-unique solutions resulting in a funny behavior. That's maybe a story for another PR or maybe just a fact of life.

[adam-urbanczyk]

Still need to add the pic, but I don't want touch the code. I think yo ucan start reviewing @jmwright .

[jmwright]

@adam-urbanczyk I think this looks good. The warning about the ftol parameter is the only thing I see that might be introducing a bug or unexpected behavior. All the rest are minor things.

[adam-urbanczyk]

OK this is it, I think. You can merge @jmwright if everything looks good.

Thanks to all for feedback and testing, thanks to @marcus7070 for the cool spindle model and to @bernhard-42 for the bearing assy that allowed me to improve the solver.

[jmwright]

Thanks for all the work you've put into this @adam-urbanczyk !

[shimwell]

Congrats on this really really nice feature