3D Printing for the Math Classroom

George Hart, April 1, 2015.  

Notes from classroom demonstrations in MfA workshop.  There are many ways to design STL files that can be built on a 3D printer. In this workshop, I’m emphasizing the RegionPlot3D function of Mathematica.  Using it, algrbraic descriptions of simple 3D regions can be created and combined with Boolean operations to create more complex forms of interest.  They are then exported to STL files for fabrication. The RegionPLot3D function is very powerful in allowing you to describe arbitrary shapes.  It is also fairly reliable for generating a buildable STL file.  (That is important because often programs for 3D design do not produce reliable STL files that can actually be built on a 3D printer.)

These notes start with the 2D analog, the RegionPlot function, before moving up to three dimensions. Execute each of the examples below (with Shift-Enter) and change them in various ways until you understand what each part of each expression does.  Also read the documentation in the Help/DocumentationCenter menu for all the functions and options used.

RegionPlot --- 2D Examples

First explore 2D examples with the RegionPlot function.  The first argument is the “characteristic function” of the region you want to show, i.e., a function you create, which is True where you want solid and False where you want air. The second and third arguments give a range for the minimum and maximum values of x and y.  Here’s a square with a hole in it.

Note:Execute each Mathematica example with Shift-Enter to see the output:

RegionPlot[x^2+y^2>1, {x,-2,2},{y,-2,2}]

Use Boolean operations to combine simple characteristic functions into more complex characteristic functions. In Mathematica, And, Or, and Not can be written as &&, ||, and !, but it is often clearer to write out Not as seen here:

RegionPlot[x^2+y^2<2^2  &&  Not[y>1], {x,-2,2},{y,-2,2}]

A trick for creating a thin region near the equation “Left=Right” is to state that the absolute value of the difference is small, i.e., Abs[Left-Right]<width.  Here is a thickened version of “y=Sin[x]”.  You can assign values to variables by using an equal sign.  Separate statements with a semicolon.

width=.2;
RegionPlot[Abs[y-Sin[x]]<width, {x,-Pi,Pi},{y,-1.5,1.5}]

If fine details of your shape aren’t being rendered, add a PlotPoints option that tells Mathematica to look harder. This slows things down if you specify very many points.  The first version below shows gaps in an intended Sine wave.  The second version fixes the gaps.  (You can type the arrow as two characters: hyphen followed by greater-than.)

width=.05;
RegionPlot[Abs[y-Sin[x]]<width, {x,-Pi,Pi},{y,-1.5,1.5}]

width=.05;
RegionPlot[Abs[y-Sin[x]]<width, {x,-Pi,Pi},{y,-1.5,1.5},PlotPoints->50]

RegionPlot3D

Everything above generalizes to 3D.  Use the RegionPlot3D function and add an argument which is the range of z.  Here’s a cube with a hole through it.  Note that you can spin the object to see all sides, by dragging it with the mouse.  Explore many examples to master this material.

RegionPlot3D[Abs[x]<1 && Abs[y]<1 && Abs[z]<1 && x^2+y^2>0.5^2,{x,-2,2},{y,-2,2},{z,-2,2}]

For greater fidelity, include the PlotPoints option, which increases the computation time: (But remember that there is no point in requesting too much fidelity, as the 3D printer has finite accuracy.)

RegionPlot3D[Abs[x]<1 && Abs[y]<1 && Abs[z]<1 && x^2+y^2>0.5^2,{x,-2,2},{y,-2,2},{z,-2,2},PlotPoints->60]

Here’s a sliced cone, revealing a hyperbola as a conic section.

RegionPlot3D[x^2+y^2<z^2  && Abs[x-.5]>0.02,{x,-1,1},{y,-1,1},{z,-1,0},PlotPoints->50]

Here’s a replacement light dimmer switch designed as a hemisphere with a 1/4 inch hole that goes in to only a certain depth.

radius=.8;
holeRadius=1/8;
depth=.6;
RegionPlot3D[x^2+y^2+z^2<radius^2  && z>0 &&  Not[x^2+y^2<holeRadius^2  &&  z<depth],{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->50]

Exporting the STL file

After executing any of the above examples, a special variable “%” holds the value of the most recently computed result.  It can be exported as an STL file with the Export command, choosing a file name (given in quotes) that ends in .STL  You may want to re-execute your command with an increased value of PlotPoints just before exporting to get a more accurate model.

Export["example.stl",%]

Alternatively, you can save the output of RegionPlot3D in a variable and export that:

result=RegionPlot3D[x^2+y^2+z^2<1^2,{x,-1,1},{y,-1,1},{z,-1,1}]

Export["ball.stl",result]

Manipulate

Mathematica has a Manipulate function which can be very handy when designing if you need to adjust a parameter and visually find a value for it that you like.  The first argument to Manipulate is the expression that includes the parameter, the second argument is the range for an automatically created slider that adjusts the parameter.  As you vary the slider, your 3D object can change.  (You can give multiple parameters and Mathematica will create multiple sliders for you.)  Slide the slider to see how it works. Here, the slider adjusts the parameter “a” which determines the eccentricity of the ellipsoid:

Manipulate[RegionPlot3D[x^2+y^2+(z/a)^2 <1 ,{x,-2,2},{y,-2,2},{z,-2,2}],{a,0.1,2}]

Egg Carton

This illustrates one way of thickening the surface given by z=f(x,y).  The points where z=Sin[x]+Sin[y] form an infinitely thin surface, so can not be built on a 3D printer.  We could build the entire solid region under the surface, but that uses a lot of plastic:

RegionPlot3D[z<Sin[x]+Sin[y],{x,-10,10},{y,-10,10},{z,-2,2},BoxRatios->Automatic, PlotPoints->70]

Instead, we can use the  Abs[difference]<thickness  trick to find the points near the surface and just build a curving shell:

thickness=0.2;
RegionPlot3D[Abs[z-(Sin[x]+Sin[y])]<thickness,{x,-10,10},{y,-10,10},{z,-2.2,2.2},BoxRatios->Automatic, PlotPoints->70]

Venn Diagram Candy Dish

We will make this twice.  The first method is ugly and inefficient, with code for the basic unit repeated three times with cut and paste.  This is motivation to learn about functions and redo the example below.  The variables below are identified here with (* comments *).  It is good programming practice to include comments which make it easy for humans to understand what you are thinking.  Text inside star parentheses is ignored, so you can add documentation.  Some blank lines can also help make complex code more readable:

r1=2;   (* outer radius *)
r2=1.6;  (* inner radius *)
top=1;   (* height at top *)
floor=0.25;  (* height of the floor *)

a1=1;b1=0;  (* xy centers of the three circles *)
a2=-1;b2=0;
a3=0;b3=Sqrt[3];

RegionPlot3D[

((x-a1)^2+(y-b1)^2<r1^2 && (x-a1)^2+(y-b1)^2>r2^2  &&  0<z<top) || ((x-a1)^2+(y-b1)^2<r1^2 && 0<z<floor) ||

((x-a2)^2+(y-b2)^2<r1^2 && (x-a2)^2+(y-b2)^2>r2^2  &&  0<z<top) || ((x-a2)^2+(y-b2)^2<r1^2 && 0<z<floor) ||

((x-a3)^2+(y-b3)^2<r1^2 && (x-a3)^2+(y-b3)^2>r2^2  &&  0<z<top) || ((x-a3)^2+(y-b3)^2<r1^2 && 0<z<floor),

{x,-3,3},{y,-2,4},{z,0,1},BoxRatios->Automatic,PlotPoints->50]

Functions

If an idea is repeated in your design, try to write a function which encapsulates that idea.  Here, there are three identical cylindrical “bowls” so let’s make a function and call it three times.  Name a fuction whatever you like, but start with a small letter so you can be sure not to redefine a built-in function.  (The built-in functions start with capital letters.)  Use square brackets, put an underbar after your dummy variables, and use colon-equals.  Here’s a function to square a variable, plus some illustrations of using the function.

f[x_]:=x^2

f[7]

f[w]

f[f[f[2]]]

You can define a function that takes two or more variables.  Separate them with commas:

add[a_,b_]:=a+b

add[3,5]

Venn Diagram Candy Dish -- Now with a Function

Take two.  Define a function for the basic bowl shape, centered at the origin.  As a characteristic function, it returns True for points inside the region we desire to build.  (The top, floor, r1 and r2 parameters were defined above.)

cylinder[x_,y_,z_]:=x^2+y^2<r1^2  &&  x^2+y^2>r2^2  &&  0<z<top ||  x^2+y^2<r1^2 && 0<z<floor

cylinder[0,0,0.1]

Now the complete candy dish can be described by caling the function three times, with translated arguments:

RegionPlot3D[cylinder[x-a1,y-b1,z]  || cylinder[x-a2,y-b2,z] ||cylinder[x-a3,y-b3,z],{x,-3,3},{y,-2,4},{z,0,1},BoxRatios->Automatic,PlotPoints->50]

Torus

Here’s a function that describes a torus.  Becaue the torus has circular symmetry, we can think about it using r and z, then later substitute r=Sqrt[x^2+y^2] to get a description in xyz coordinates. In r and z, we just need to say we are interested in points inside the circle of radius r2 centered at (r1,0), so (r-r1)^2 + z^2 < r2^2.  Substituting to eliminate r gives:

r1=4;
r2=1;
torus[x_,y_,z_]:=(Sqrt[x^2+y^2]-r1)^2+z^2<r2^2

We can check this gives a single torus :

RegionPlot3D[torus[x,y,z], {x,-8,8},{y,-8,8},{z,-4,4},BoxRatios->Automatic]

Links of a Chain

Let’s use the avove function to put two tori together in a linked configuration.  We can rotate one 90 degrees about the x axis by swapping y and z.  The other is translated in the x direction by subtracting 3 from x:

RegionPlot3D[torus[x,z,y]||torus[x-3,y,z], {x,-8,8},{y,-8,8},{z,-5,5},BoxRatios->Automatic]

Polyhedra

Flat surfaces of polyhedra can be defined by taking the dot product of the vector {x,y,z} with a vector in the normal direction.  Here’s a regular octahedron, expressed using a dot product.  The vector {x,y,z} is dotted with a direction vector, {a,b,c} to give a measure of how far {x,y,z} projects in that direction.  Use the directions to the corners of a cube to make its dual, the octahedron.

d=1;
RegionPlot3D[ 
   -d<{x,y,z}.{1,1,1}<d && 
-d<{x,y,z}.{-1,1,1}<d&& 
-d<{x,y,z}.{1,-1,1}<d&& 
-d<{x,y,z}.{-1,-1,1}<d,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->50]

To better understand how it works, try changing d, e.g., the same code but with d=2 gives a cuboctahedron. (Also try roughly 2.4, for regular octagons.)

d=2;
RegionPlot3D[ 
   -d<{x,y,z}.{1,1,1}<d && 
-d<{x,y,z}.{-1,1,1}<d&& 
-d<{x,y,z}.{1,-1,1}<d&& 
-d<{x,y,z}.{-1,-1,1}<d,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->50]

The above generalizes to any set of directions, e.g., you could make a dodecahedron from its normal directions (which are easy to work out---good exercise).  For the octahedron, because it has the same symmetry as the xyz axes, we can do a shortcut and distribute the Absolute value function over the components of {x,y,z}:

d=1;
RegionPlot3D[ 
   Abs[{x,y,z}].{1,1,1}<d ,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->50]

A

Using a Table to Create Multiple Objects

If you want to produce a series of objects with one command, the Table[] function is very useful.  Its second argument is an iterator that specifies te range of a variable.  The first argument to Table is a function of the variable that you want to evaluate for each value of the variable.  Table creates a list of all the results.  In first example below, the function us just the variable i, so we see the eleven values it takes on from 0 to 10.  In the second example, we compute i^2 for these values.  In the third example, we create eleven characteristic functions; each is for a sphere with a different z value for its center.

Table[i,{i,0,10}]

Table[i^2,{i,0,10}]

Table[x^2+y^2+(z-i)^2<1^2,{i,0,10}]

Suppose we want to create all eleven spheres described by the list of equations above.  We need to make one characteristic function which is the Or of them all.  The Apply function take a function and a list as arguments and applies the function to the list.  For example we could Apply Plus to a list of thinge we want to add and the result is an expression for their sum:

Apply[Plus,{x,y,1}]

We will Apply Or to our list of individual characteristic functions to create a single large equation for the union of the spheres:

Apply[Or,Table[x^2+y^2+(z-i)^2<1^2,{i,0,10}]]

Now we can use that characteristic function in RegionPlot3D to create all the spheres.  You can change n to change the number of spheres.

n=10;
RegionPlot3D[Apply[Or,Table[x^2+y^2+(z-i)^2<1^2,{i,0,n}]],
{x,-1,1},{y,-1,1},{z,-1,n+1},BoxRatios->Automatic,PlotPoints->40]

As another example, let’s arrange the spheres in a circle.  We’ll create a variable theta to be n equally spaced angles around a circle.  We need theta to step by 2Pi/n, so we specify that in the iterator.  The values r Cos[theta] and r Sin[theta] are the centers of the n spheres.

n=12;
r=6;RegionPlot3D[Apply[Or,Table[(x-r Cos[theta])^2+(y-r Sin[theta])^2+z^2<1^2,{theta,2Pi/n,2Pi,2Pi/n}]],
{x,-10,10},{y,-10,10},{z,-1,1},BoxRatios->Automatic,PlotPoints->40]

Rotations

We saw how to rotate 90 degrees (by swapping two of x, y, and z) in the torus function example above, but for an arbitrary rotation, it is handy to use a rotation matrix.  A built-in function can create a 3-by-3 rotation matrix for any given axis and angle to rotate about that axis.  Here’s one for rotating 45 degrees about the X axis:

m=RotationMatrix[Pi/4,{1,0,0}]

Multiply the matrix times the vector {x,y,z} with the dot to indicate matric-vector multiplication.  This gives a rotated vector:

m.{x,y,z}

We' ll try it out by creating a cylinder function which takes a vector as its argument.  First plot the cylinder in its standard orientation.  Then with “m.” we create a rotated vector of arguments (so the cuylinder we see is rotated 45 degrees the opposite way.)

cylinder[{x_,y_,z_}]:=x^2+y^2<1 && Abs[z]<.25

RegionPlot3D[cylinder[{x,y,z}],{x,-2,2},{y,-2,2},{z,-2,2},PlotPoints->30]

RegionPlot3D[cylinder[m.{x,y,z}],{x,-2,2},{y,-2,2},{z,-2,2},PlotPoints->30]

Here' s an example combining three different rotations in a Table, to make three superposed disks rotated 0, 45 and 90 degrees:

RegionPlot3D[Apply[Or,Table[cylinder[RotationMatrix[th,{1,0,0}].{x,y,z}],{th,0,Pi/2,Pi/4}]],{x,-2,2},{y,-2,2},{z,-2,2},PlotPoints->30]

Tubes

There are other 3D functions in Mathematica that might be useful in special situations, but be warned that some are buggy and may produce an STL file that can not be 3D printed.  If your desired shape is a tube, like a torus, a knot, or a 3D Lissajous figure, try ParametricPlot3D.  The Tube[] PlotStyle lets you specify a desired radius for the tube.  Here’s a random 3D Lissajous figure.

ParametricPlot3D[{Cos[5t],Sin[3t],Cos[2t]},{t,0,2 Pi},PlotStyle->Tube[.1],PlotRange->All]

Note that as of this writing, there is a bug in Mathematica which puts caps on that face the wrong direction.  Thic can be detected with a program like VisCam View. It can be repaired with a program like NetFabb.

Graphics3D

Here’s another approach to generating 3D consructions, not using RegionPlot3D.  The Graphics3D function takes a list of predefined shapes, including sphere, cuboid, cylinder.  See the documentation for a description of their arguments.  Here are a few examples to illustrate the idea.  But be aware that when you export these things to an stl file, interior surfaces will be present where obhects overlap.  This violates the official rules of legal stl files and some 3D printers may have trouble with the file.

Graphics3D[{Sphere[{0,0,1},1],Sphere[{0,3,0},.5],Cuboid[]}]

Graphics3D[{Table[Sphere[{0,0,i},1],{i,0,5}], Cylinder[{{0,0,0},4{1,1,1}},.5]}]

More Examples

Some examples you can practice with, which I didn’t discuss in class:

Here is a “superellipse” as defined by the Danish designer Piet Hein.  Its equation is analogous to that of a sphere, but changing the exponent from 2 to 2.8.  (You can buy metal superellipse “ice cubes” of this shape that you store in the freezer then use to keep drinks cold without diluting them.)  Larger exponents make an object closer to a cube.  Also try it to see what results when the exponent is in the range 0.5 to 1.5.

c=2.8;
RegionPlot3D[Abs[x]^c+Abs[y]^c+Abs[z]^c<1,{x,-1,1},{y,-1,1},{z,-1,1}]

Here's a visual 10-by-10 multiplication table.  The height of each bar corresponds to the product of the x and y positions.

RegionPlot3D[z<Floor[x]*Floor[y],{x,1,11},{y,1,11},{z,0,100},PlotPoints->100]

Compare it to a visual 10-by-10 addition table, where the height of each bar corresponds to the sum of the x and y positions.  Here, every staircase has the same slope.

RegionPlot3D[z<Floor[x]+Floor[y],{x,1,11},{y,1,11},{z,0,20},PlotPoints->100]

Try making a nut and a bolt that can screw together.  The helix here is based on thinking of z a free parameter, choosing the center of the helix as x=r Cos[a z] and y=Sin[a z] for each z, then selecting points not further than epsilon from that center.  The bolt is the helix plus a cylindrical core.

r=.2;
epsilon=.2;
a=10;
RegionPlot3D[(x-r Cos[a z])^2+(y-r Sin[a z])^2<epsilon^2  || x^2+y^2<r^2,{x,-.5,.5},{y,-.5,.5},{z,-1.5,1.5},PlotPoints->75,BoxRatios->Automatic]

The nut is the complement of the bolt, so just put a Not[] around everythig and shorten it.

r=.2;
epsilon=.2;
a=10;
RegionPlot3D[Not[(x-r Cos[a z])^2+(y-r Sin[a z])^2<epsilon^2  || x^2+y^2<r^2],{x,-.5,.5},{y,-.5,.5},{z,0,.5},PlotPoints->50,BoxRatios->Automatic]

Animations

In addition to stl files, you might want to export images of your designs of animations.  Just use a file name with any standard three letter extension and Mathematica will convert to that format, e.g., to make an svg image file of something you just created, use:

Export["my-thing.svg",%]

If you want to make an animation, make a table of images by having an expression with a parameter as the first argument to Table and a range with an increment as the second argument.  In this example, the parameter c will take on values 0.0, 0.1, 0.2, ... 0.9, 1.  Then export the table as an AVI file for video or an animated GIF file for a web page.  (I’m creating the table just once and saving in in a variable so I can export it twice.)

table=Table[RegionPlot3D[x^2+y^2<z^2  &&  x<c ,{x,-1,1},{y,-1,1},{z,-1,0}],{c,0,1,.1}];

Export["anim.avi",table]

Export["anim.gif",table]