Coding 3D parts in Python using FreeCAD

For quite a long time I have been a happy user of OpenSCAD software, mostly for designing 3D printed parts. For a recent project, I needed to create a complex geometry, and my starting point was a data file, so it made a lot of sense to that with a program. 

Having been using OpenSCAD for a long time, I got the job done with it but creating the STL file took more than 10 hours. I could live with that, but I was not looking for an STL file as an output but a STEP or BREP file instead. So I started reviewing what I could do in Python using FreeCAD. It turned out that I could do a lot, but it was not as easy as using OpenSCAD. And it was not very fast either, but I could get it working in the end.

The code below is what I used to create the model of the image of the cube, but for reasons unknown, all the shapes are hidden (not visible) when I open the file with FreeSCAD.

It is very simple but, although there is a lot of documentation out there, I did not find all of it when I was looking for examples. That is why I thought it might be of use for other users.

import FreeCAD
import Part
from FreeCAD import Base, Vector

lista = [[(0,0,0),(0,50,0)], \
        [(0,50,0),(50,50,0)], \
        [(50,50,0),(50,0,0)], \
        [(50,0,0),(0,0,0)], \
        [(0,0,0),(0,0,50)], \
        [(0,50,0),(0,50,50)], \
        [(50,0,0),(50,0,50)], \
        [(50,50,0),(50,50,50)], \
        [(0,0,50),(0,50,50)], \
        [(0,50,50),(50,50,50)], \
        [(50,50,50),(50,0,50)], \
        [(50,0,50),(0,0,50)], \
        ]

def crossproduct(first, other=FreeCAD.Vector(0,0,1)):
    """crossproduct(Vector,Vector) - returns the cross product of both vectors. 
    If only one is specified, cross product is made with vertical axis, thus returning its perpendicular in XY plane"""
    if isinstance(first,FreeCAD.Vector) and isinstance(other,FreeCAD.Vector):
        return FreeCAD.Vector(first.y*other.z - first.z*other.y, first.z*other.x - first.x*other.z, first.x*other.y - first.y*other.x)

DOC = FreeCAD.newDocument()

i=1
points=set()
for a,b in lista: 
    points.add(tuple(a))
    points.add(tuple(b))
    b = Vector(b)
    a = Vector(a)
    c = b - a
    o2 = DOC.addObject('Part::Cylinder','cilindro'+str(i))
    o2.Height = c.Length
    if c.x!=0 or c.y!=0:
        c.normalize()
        u=crossproduct(c) #with (0,0,1)
        u.normalize()
        v=crossproduct(c,u)
        v.normalize()
        m = App.Matrix( u[0],v[0],c[0],a[0], u[1],v[1],c[1],a[1], u[2],v[2],c[2],a[2], 0,0,0,1)
        o2.Placement=FreeCAD.Placement(m)
    else:
        o2.Placement=FreeCAD.Placement(Vector(a),FreeCAD.Rotation(0,0,0),Vector(0,0,0))
    i+=1

i=1
for x,y,z in points:
    o = DOC.addObject('Part::Sphere','esfera'+str(i))
    o.Radius = 3
    o.Placement=FreeCAD.Placement(Vector(x,y,z),FreeCAD.Rotation(0,0,0),Vector(0,0,0))
    i+=1


DOC.recompute()
DOC.saveAs(u"your_own_path/cube.FCStd")