Goldberg polyhedron macro failure

Need help, or want to share a macro? Post here!
Forum rules
Be nice to others! Respect the FreeCAD code of conduct!
wangnick
Posts: 19
Joined: Thu Nov 03, 2022 9:27 am

Goldberg polyhedron macro failure

Post by wangnick »

Dear all,

I'm trying to create a Goldberg polyhedron using a Python macro in FreeCAD. Attached my attempts.

What I recognise is that the pentagon faces are fine but the hexagon faces that I'm trying to create are not planar. For a Goldberg with radius 5mm and division level 2 (geodesic(5,2,True)) I observe that some of the three vertices of some of the hexagons have a distance to the plane created by the first three vertices of about 0.05mm.

I did succeed earlier to create a proper solid by using Part.makeFilledFace(). However, that seemingly creates a curved face, which makes FreeCAD very very slow with higher-division Goldbergs (for my current design project I need a Goldberg with division level 5).

Does anyone have an idea how to adapt my macro such that the hexagons become planar? I would not care if some of the vertices are not exactly located on the circumsphere. But I'm really stretching my geometric math capabilities here ...

Kind regards,
Sebastian

PS: FYI, the core of the macro, the creation of the geodesic sphere, is a slightly adapted version of the code from the "Pyramids and Polyhedrons" workbench of Eddy Verlinden.
Attachments
goldbergtest03.FCMacro
(6.61 KiB) Downloaded 22 times
goldbergtest03.FCStd
(35.23 KiB) Downloaded 18 times
goldbergtest03.png
goldbergtest03.png (258.19 KiB) Viewed 1221 times
edwilliams16
Veteran
Posts: 3117
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Goldberg polyhedron macro failure

Post by edwilliams16 »

Maybe https://forum.freecadweb.org/viewtopic. ... 57#p214157 would be helpful.

If some of the points are 0.05 mm out-of plane, that is huge for a rounding error!
wangnick
Posts: 19
Joined: Thu Nov 03, 2022 9:27 am

Re: Goldberg polyhedron macro failure

Post by wangnick »

edwilliams16 wrote: Fri Nov 18, 2022 6:43 pm If some of the points are 0.05 mm out-of plane, that is huge for a rounding error!
No, I don't think that is a rounding error, I think that is a problem in my construction algorithm. But I am not sufficient proficient in polyhedric geometry to correct that problem. Any help appreciated.

What's worse, the Cura slicer is not throwing an error when trying to slice my .stl file, and that might well be caused by my very slightly curved faces:

Code: Select all

2022-11-18 19:14:43,140 - DEBUG - [EngineOutputThread] UM.Backend.Backend._backendLog [115]: [Backend] [2022-11-18 19:14:43.142] [warning] Mesh has overlapping faces!
...
2022-11-18 19:14:44,704 - DEBUG - [EngineOutputThread] UM.Backend.Backend._backendLog [115]: [Backend] [2022-11-18 19:14:44.708] [warning] Failing to discretize parabola! Must add an apex or one of the endpoints.
2022-11-18 19:14:46,302 - ERROR - [MainThread] CuraEngineBackend.CuraEngineBackend._onBackendQuit [964]: Backend exited abnormally with return code 3221225477!
Kind regards,
Sebastian
edwilliams16
Veteran
Posts: 3117
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Goldberg polyhedron macro failure

Post by edwilliams16 »

The geodesic sphere is presumably good, so the problem is taking the dual. Where are you getting the dual algorithm from?
edwilliams16
Veteran
Posts: 3117
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Goldberg polyhedron macro failure

Post by edwilliams16 »

I didn't try to debug your dual, but mine seems to work. Your dual vertices weren't correct. I found the dual vertices, then hacked the edges, wires, faces somewhat inefficiently

Code: Select all

import FreeCAD as App
from FreeCAD import Vector as V3

doc = App.ActiveDocument

def dualofsphtriangle(rad, n1, n2, n3):
    ''' Dual vertex location from triangle'''
    return r *(n1.cross(n2) + n2.cross(n3) + n3.cross(n1))/n1.dot(n2.cross(n3))

shp = doc.getObject("GoldbergGeodesicSphere").Shape
r = shp.Vertexes[0].Point.Length # 5 mm

objs = []
dpts = []
for w in shp.Wires:
    pts = [v.Point.normalize() for v in w.Vertexes]
    dpt = dualofsphtriangle(r, *pts)
    dpts.append(dpt)
    obj = Part.show(Part.Vertex(dpt))
    obj.ViewObject.PointSize = 4
    objs.append(obj)

cmpd = doc.addObject("Part::Compound","DualVertices")
cmpd.Links = objs

doc.recompute()

def wireOnDualPlane(v, dpts ):
    ''' find wire dual to vertex v '''
    n = v/v.Length
    pts = [pt for pt in dpts if (v - n.dot(pt) * n).Length < 1e-7]
    tempEdges = []
    for i, v1 in enumerate(pts):
        for j, v2 in enumerate(pts):
            if j > i:
                tempEdges.append(Part.LineSegment(v1, v2).toShape())

    sortEdges = sorted(tempEdges, key = lambda e: e.Length)
    #either 5 or 6 shortest edges - definitely a hack
    if abs(sortEdges[4].Length - sortEdges[5].Length) < 1e-7:
        culledEdges = sortEdges[0:6]
    else:
        culledEdges = sortEdges[0:5]
    return Part.Wire(*Part.sortEdges(culledEdges))


faces = []
for v in shp.Vertexes:
	faces.append(Part.Face(wireOnDualPlane(v.Point , dpts )))

shell= Part.Shell(faces)
Part.show(shell, "GoldbergPolyhedronShell")

Screen Shot 2022-11-18 at 5.17.51 PM.png
Screen Shot 2022-11-18 at 5.17.51 PM.png (20.89 KiB) Viewed 1065 times
Attachments
goldbergtest03_ew.FCStd
(110.06 KiB) Downloaded 24 times
wangnick
Posts: 19
Joined: Thu Nov 03, 2022 9:27 am

Re: Goldberg polyhedron macro failure

Post by wangnick »

edwilliams16 wrote: Fri Nov 18, 2022 11:56 pm The geodesic sphere is presumably good, so the problem is taking the dual. Where are you getting the dual algorithm from?
The algorithm is of my doing, inspired by the Polar Reciprocation section of https://en.wikipedia.org/wiki/Dual_polyhedron, and by the Blender add_mesh_geodesic_domes source code (class mesh method dual in vefm_271.py). All mistakes are mine.
edwilliams16 wrote: Sat Nov 19, 2022 3:28 am Your dual vertices weren't correct.
I'll have a look. Thanks for investigating!

Kind regards,
Sebastian
wangnick
Posts: 19
Joined: Thu Nov 03, 2022 9:27 am

Re: Goldberg polyhedron macro failure

Post by wangnick »

edwilliams16 wrote: Sat Nov 19, 2022 3:28 am I found the dual vertices ...
Ed, I see you used the formula

Code: Select all

(n1.cross(n2) + n2.cross(n3) + n3.cross(n1))/n1.dot(n2.cross(n3))
and do some magic using the normalised vertex vectors of the triangular face of the geodesic sphere. I haven't seen that before. Which kind of point does that compute? What is the geometric reasoning behind that? Super cool, and I'd like to understand how that assures coplanarity of the hexagons ...

Kind regards,
Sebastian
edwilliams16
Veteran
Posts: 3117
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Goldberg polyhedron macro failure

Post by edwilliams16 »

wangnick wrote: Sat Nov 19, 2022 2:41 pm
The algorithm is of my doing, inspired by the Polar Reciprocation section of https://en.wikipedia.org/wiki/Dual_polyhedron, and by the Blender add_mesh_geodesic_domes source code (class mesh method dual in vefm_271.py). All mistakes are mine.
If you used the Dorman-Luke construction of the above, it won't work because the Goldberg Polyhedron isn't uniform - meaning the vertices are not all equivalent - in that you can't rotate one into any other while preserving the shape. I had to go back to the more general polar reciprocation definition. I couldn't find a reference that spelled it out using vectors in 3d (which is why I asked you), so I gave up and derived it. I'll write it up if you like. One thought I had was that a triangular face numerically unambiguously maps into a dual vertex. Faces with more edges (like you would face retrieving the Geodesic Sphere as the dual of the Goldberg polyhedron) provide an over-determined system of equations for its dual vertex. But, even with triangular faces there is still floating point error, but FreeCAD's tolerances seems to accommodate those OK, so maybe that isn't a problem. Just pick three (not collinear) vertices of the face and go with them.
What would be a real improvement algorithmically would be to transfer the topology of the polytope to its dual more efficiently - ie connecting up the dual vertices with edges in a less of a brute force algorithm. Right now I loop through the vertices of the polytope and find the dual vertices of the entire polytope that lie on the plane dual to the looped vertex. These are the vertices of the face dual to the looped vertex. Then I have to order them to make a wire. The method I used I suspect wouldn't generalize to less regular polytopes.
edwilliams16
Veteran
Posts: 3117
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: Goldberg polyhedron macro failure

Post by edwilliams16 »

Here's the short version.

The plane dual to a vertex at p w. r. t. a sphere of radius r is given by the set {x| p.x = r^2}

If we have a triangular face with vertices p1, p2 and p3, then the vertex dual to the face is at the intersection of the three planes dual to the three vertices. So the vertex dual to the face, x, satisfies:
eq1.png
eq1.png (9.88 KiB) Viewed 877 times
You can verify by inspection the solution for x is (Cramer's rule in vector notation) given by
eq2.png
eq2.png (10.77 KiB) Viewed 877 times
which give or take normalizations is the equation you were wondering about. I'm sure it has been known for a century, but I failed to google it.
I used the circumsphere of the Geodesic sphere as the reference sphere, so the planes dual to the vertices pass through the vertices, normal to the radius.
wangnick
Posts: 19
Joined: Thu Nov 03, 2022 9:27 am

Re: Goldberg polyhedron macro failure

Post by wangnick »

Superb, Ed. I feel I start to understand Goldberg polyhedra a bit better.

As to the transfer the topology to its dual I tried to improve your macro a bit by extracting and then using something of a data structure:

Code: Select all

import FreeCAD as App

doc = App.ActiveDocument

def dualofsphtriangle(n1, n2, n3):
    ''' Dual vertex location from triangle'''
    return (n1.cross(n2) + n2.cross(n3) + n3.cross(n1))/n1.dot(n2.cross(n3))

shp = doc.getObject("GoldbergGeodesicSphere").Shape
r = shp.Vertexes[0].Point.Length # 5 mm

#objs = []
dpts = dict()
vd = dict() # dict mapping vertex to set of attached wires. vertex is tuple of cartesian coordinates. wire is frozenset of vertexes (again as tuples of cartesian coordinates).
for w in shp.Wires:
    we = frozenset(tuple(round(xyz,9) for xyz in v.Point) for v in w.Vertexes)
    for ve in we: vd.setdefault(ve,set()).add(we)
    pts = [v.Point.normalize() for v in w.Vertexes]
    dpt = dualofsphtriangle(*pts)
    dpts[we] = dpt

dvmax = max(dpt.Length for dpt in dpts.values())
for dpt in dpts.values(): 
    dpt.multiply(r/dvmax)
    #obj = Part.show(Part.Vertex(dpt))
    #obj.ViewObject.PointSize = 4
    #objs.append(obj)

#cmpd = doc.addObject("Part::Compound","DualVertices")
#cmpd.Links = objs

doc.recompute()

faces = []
for ve,wes in vd.items():
    # We build a dual face for each vertex, using the dual vertex positions of the wires attached to the vertex. 
    # For each vertex we collected a set of wires above which we now bring into a proper sequence.
    pwe = wes.pop()
    lpts = [dpts[pwe]]
    while wes:
        ves = set(pwe) # The vertexes of the previous wire (each a tuple of cartesian coordinates) as a modifyable set
        ves.remove(ve) # We want to find a wire we in the remaining set of wires that shares a second vertex with the previous wire
        we = None
        for wei in wes:
            if any(ve in ves for ve in wei):
                we = wei
                break
        if we is None: break
        wes.remove(we)
        lpts.append(dpts[we])
        pwe = we
    if wes:
        print("Error: Failed to bring wires attached to vertex %s in sequence"%str(ve))
    else:
        faces.append(Part.Face(Part.makePolygon(lpts+lpts[:1])))

shell= Part.Shell(faces)
Part.show(shell, "GoldbergPolyhedronShell")
It does not use edge information and so investigates a few vertices too many but seems to do the job.

Kind regards,
Sebastian
Post Reply