A friend wanted to 3D-print a shape to demonstrate using calculus to find the volume of solids of known cross-section.
The shape he wanted was a graph of $sin(x)$ vs $x^2$,
where each vertical slice of the intersection was a square.
Here’s the graph, with $y_1 = sin(x)$ and $y_2 = x^2$. The blue lines show the edge of each square.
He couldn’t figure out how to do this in a CAD program (I’m not even sure if it’s possible),
so he asked me if I could write some code to render it.
This post was written in literate Python. Download
buildshape.py and run it yourself!
The first attempt I made was written directly in OpenSCAD;
I wrote a small loop to union together a collection of cube()
calls, one for each slice.
This had several problems, however.
First, because I was using cubes it looked very chunky unless the slices were very small.
Second, I ran into a bug in OpenSCAD
which causes something like $O(n^2)$ performance for this scenario
and the file took about 2 hours to render.
Finally, the STL file it emitted was 1.8 MB and crashed the slicer.
Clearly, this solution was not going to work.
Instead, I wrote a python script to output an OpenSCAD file containing a single polyhedron,
which I could then render into an STL file and hand off to him to slice and print.
In addition to not crashing the slicer,
this approach also had the advantage of allowing the resolution to be much coarser
while still avoiding the ‘stair-step’ problem of the original multiple cube()
approach.
Setting up
First, the formulae:
40 41
| def y1(x): return math.sin(x) def y2(x): return x*x
|
$sin(x) = x^2$ is true for $x = 0$ and $x \approx 0.876$, so that’s where I’ll start and stop.
44 45
| START = 0 END = 0.876726
|
I’ll draw 50 slices, and scale everything up by 30.
48 49
| SLICES = 50 SCALE = 30
|
OpenSCAD’s polyhedron()
function
takes two arguments, points
and triangles
.
(Versions 2014.03 and later can also take faces
instead of triangles
, but I’m still using 2014.01.)
points
is a list of $[x, y, z]$ triplets, and triangles
then indexes into the list of points to build the triangles.
I therefore make a function point(x, y, z)
which will save the point into the list and return the index into this list.
Most of the program then works using these indexes, rather than the actual point coordinates.
58 59 60 61
| points = [] def point(x, y, z): points.append( (x*SCALE,y*SCALE,z*SCALE) ) return len(points)-1
|
Slicing the shape

The first step is to slice the shape and get a list of squares.
(OpenSCAD expects the points on a triangle to be clockwise when looking at them from the outside,
so throughout the program I maintain this clockwise orientation.) Here's the square; line $\overline{P_2 P_3}$ is the part of the square that sits on the graph. The remainder of the square sits above it.
I slice the shape into
SLICES
slices, and output the above square for each slice.
This is also where I convert the coordinates into point-list indexes as explained above.
70 71 72 73 74 75 76 77
| def get_slices(): for x in scipy.arange(START, END, (END-START)/SLICES): yield ( point(x, y1(x), y1(x)-y2(x)), point(x, y2(x), y1(x)-y2(x)), point(x, y2(x), 0), point(x, y1(x), 0), )
|
Building the polyhedron
All of these slices I have are inside the shape,
which isn’t the part that’s visible.
The next step is to get all of the faces of the polyhedron
by connecting the edge of adjoining slices into quadrilaterals.
First, I make an “end cap” using the slice at the start of the shape.
This is necessary because otherwise the end will be left ‘open’ if the first slice is not of size 0,
and the shape cannot be exported
because it is not a valid 2-manifold.
Next I generate quadrilaterals conecting each of the sides of ajoining slices.
Doing this is relatively simple:
$$
\mbox{Given slices}
\substack {
\square T_0 T_1 T_2 T_3 \\
\square P_0 P_1 P_2 P_3
}
\mbox{, emit faces}
\substack {
\square T_0 P_0 P_1 T_1 \\
\square T_1 P_1 P_2 T_2 \\
\square T_2 P_2 P_3 T_3 \\
\square T_3 P_3 P_0 T_0
}.
$$
(Notice that they all go clockwise when facing outward.)
113 114 115 116 117 118 119 120 121 122
| for this in slices[1:]: def face(a, b): return ( this[a], prev[a], prev[b], this[b], ) yield face(0, 1) yield face(1, 2) yield face(2, 3) yield face(3, 0) prev = this
|
The last step in making the polyhedron is to cap off the back end the same as I did with the front end (see above).
Convert quadrilaterals into triangles
I now have a list of faces for the polyhedron. For newer versions of OpenSCAD I could just pass this list in to the polyhedron()
function; unfortunately as I mentioned before my version of OpenSCAD only accepts a list of triangles.
Therefore, I take each face (a quadrilateral $\square F_0 F_1 F_2 F_3$) and divide it into two triangles, $\triangle F_0 F_1 F_3$ and $\triangle F_1 F_2 F_3$. (Again, I’m being very careful to keep them in clockwise order.)
130 131 132 133
| def get_triangles(): for face in get_faces(): yield (face[0], face[1], face[3]) yield (face[1], face[2], face[3])
|
Write the file
Finally, I print out the OpenSCAD file, which simply consists of a single call to polyhedron()
:
137 138 139 140 141 142 143 144 145 146 147
| print "polyhedron(" print " convexity = 1," print " triangles = [" for triangle in get_triangles(): print " [{0}, {1}, {2}],".format(*triangle) print " ]," print " points = [" for point in points: print " [{0}, {1}, {2}],".format(*point) print " ]" print ");"
|
And finally, here’s what the shape looks like when it’s printed out.
In addition to the Python program linked at the top of this post,
you can also download
the OpenSCAD file and
the STL file.
