Mathematical landmines

[neri-engineering]

I realize that these are OpenCascade problems and that the math must be very complicated, but nonetheless I keep finding these.  I must dance around issues very frequently.  Here is short code to generate the explosive in question:

import cadquery as cq

globe = cq.Workplane("XY")

globe = globe.sphere(10.125)

teh_box = cq.Workplane("XZ")

teh_box = teh_box.box(13.982755172285637, 21.2625, 21.2625, centered=True)

globe = globe.intersect(teh_box)

globe = globe.edges().fillet(2.4494254642267013)

# Image 1 shows 'globe'.

bore = cq.Workplane("XY")

bore = bore.circle(9.426271240411136)

bore = bore.extrude(21.2625)

bore = bore.translate((0, 0, -10.63125))

bore = bore.intersect(teh_box)

bore = bore.edges("|Z").fillet(2.7442687282902405)

# Image 2 shows 'bore'.

part = globe.intersect(bore)

# Image 3 shows 'part'.

cyl = cq.Workplane("XY")

cyl = cyl.circle(5.0)

cyl = cyl.extrude(until= 16, both= True)

garbage = part.union(cyl)

# Image 4 shows 'garbage'.

show_object(garbage)

The four images alluded to in the code are as follows, in order:

Now on to another insight that I've had.  It's regarding orthogonal rotations, namely rotations along X, Y, or Z axes, where the angle of rotation of some multiple of 90°.

The matrices which are created as a result oftentimes have terms such as 'sin(pi/2)' and 'cos(pi/2)'.  Unfortunately these come close to the desired value but because we cannot express π/2 exactly, we oftentimes get terms which ARE NOT 0 and/or 1, but very close to it.

In other words rotating your object by 90° may or may not cause axis-alignment after the fact.  I see this as yet more landmines to dance around.

There is a reason why "degrees" were chosen as units to specify rotations at the high level API.  So good for CadQuery on making this correct choice, in my opinion.  The reason why it's a good choice is because 180 is divisible by two, and 90 is divisible by 3 and by 2.  We can do whole number math for the most important rotations.

We have other interesting phenomena besides the following which I've already mentioned:

cos(90°) = 0

sin(90°) = 1

Namely we have:

cos(45°) = sin(45°)  which can be expressed with a sqrt(2) value somewhere in the formula.

But we also have:

sin(30°) = 0.5

cos(60°) = 0.5

At the very least, the orthogonal rotations having multiples of 90° should return perfect results, but after looking at the code, it seems that more landmines result.

To avoid these landmines special-case code such as this could be employed; I've now resorted to including this in all of my code.

from math     import floor      as floor

from cadquery import Assembly   as Assembly

from cadquery import Shape      as Shape

from cadquery import Workplane  as Workplane

from OCP.gp   import gp_Trsf

# Returns an angle equivalent to the input angle but in the range [0,360).

# Negative and positive input values accepted.  E.g., -21.5 -> 338.5.

def deg_canonical_360(deg_theta):

    return deg_theta - floor(deg_theta / 360) * 360

# Convenience function for rotating a piece around the x axis, right hand rule.

# The piece can be a Workplane, a Shape, or an Assembly, in CadQuery.  For

# Workplane and Shape, orthogonal rotations are handled as special cases in

# order to provide mathematical exactness.

def rotateX(piece, degrees):

    if isinstance(piece, Assembly):

        raise Exception("Assembly not implemented yet")

    canonical = deg_canonical_360(degrees)

    xform = None

    if canonical ==   0:  xform = gp_Trsf()  # Identity transform.

    if canonical ==  90:  xform = _get_x_rot_90()

    if canonical == 180:  xform = _get_x_rot_180()

    if canonical == 270:  xform = _get_x_rot_270()

    if xform is not None:

        print("SPECIAL CODE TRIGGERED")

        if isinstance(piece, Workplane):

            return piece.newObject(

                [

                    obj._apply_transform(xform)

                    if isinstance(obj, Shape)

                    else obj

                    for obj in piece.objects

                ]

            )

        if isinstance(piece, Shape):

            return piece._apply_transform(xform)

    return piece.rotate(axisStartPoint=(0,0,0),

                        axisEndPoint=(1,0,0),

                        angleDegrees=degrees)

def _get_x_rot_90():

    xform = gp_Trsf()

    xform.SetValues( 1,  0,  0,  0,

                     0,  0, -1,  0,

                     0,  1,  0,  0 )

    return xform

# Draw an 'L' as viewed from -y direction.

piece1 = Workplane("XY")

piece1 = piece1.box(0.1, 0.1, 4, centered= False)

piece2 = Workplane("XY")

piece2 = piece2.box(2, 0.1, 0.1, centered= False)

pieces = piece1.union(piece2)

pieces = rotateX(pieces, 90)

show_object(pieces)

Cheers.  Nerius.

Sent with Proton Mail secure email.