3D rigid body physics engine

Is there a package in Mathematica that will allow me to simulate collisions of 3d rigid bodies?

If not, what known libraries could I use and how?

For example, one problem I want to simulate is a coin toss. I would need to specify an infinite plane, a cylinder with an initial linear and angular momentum, and properties like friction and restitution. After running the simulation I would want to check (no need for graphics here) if it landed head or tails (or side).

As I said, I don't need the simulation to show graphics.

Edit: Someone mentioned UnityLink. Any resources about getting started will be appreciated!

• We'll need a bit more info... Are you interested in collisions? In forces between the bodies? Vibrations? What kind of problem are you trying to solve? Apr 17, 2019 at 4:26
• @b3m2a1 For now, I'm interested in collisions. Also I'm interested in simulating cubes, spheres, cilynders.
– Ivan
Apr 17, 2019 at 4:30
• There is no built-in physics engine, unfortunately. Apr 17, 2019 at 4:33
• @C.E. does Unity have such a thing? There's now UnityLink which would make it easy to work with theirs. Apr 17, 2019 at 4:37
• @C.E. True but if the OP wants to invest the time I bet it’d pay off. Apr 17, 2019 at 6:56

Since this answer was originally posted, the Blender software package has undergone significant revision causing the Python scripts to break. I have modified the scripts so that they should work from v2.79b to v2.93LTS.

Updated Blender scripts

Note: I added an update below to import position and orientation transformations and view the 3D simulation results in Mathematica.

I have used the free program Blender v2.79b to simulate the handling of 100s of complex shapes through a geometrically complex industrial machine with many moving parts including vibrating elements. So, it should be able to handle a "coin flip". I believe Blender still uses the Bullet Physics Engine as its solver. I should warn you that collision simulation can get difficult and that there are a lot of tricks of the trade you must learn to be accurate and fast in a general case.

Blender has a python interface and it can be run as a background task (Bullet also has a python interface, but I am not familiar with its operation). Since Mathematica can create text files with StringTemplate and can execute system commands, we should be able to create a python script to drive a Blender simulation.

Python Script Generation

Blender has a fairly well-documented API and there many resources one can find online to generate python script.

import bpy
from math import pi

version = bpy.app.version

for o in bpy.data.objects:
if version < (2, 80, 0):
if o.type == 'MESH' or o.type == 'EMPTY':
o.select = True
else:
o.select = False
else:
o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

if version < (2, 80, 0):
else:
bpy.ops.transform.resize(value=(1, 1, 0.1))
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "CYLINDER"
bpy.context.object.rigid_body.friction = 0.25
bpy.context.object.rigid_body.restitution = 0.75
boxObj.name = "Coin"
# Set reference to the coin
coin = bpy.data.objects["Coin"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

sce.frame_set(4)
# Translate up a little
coin.location.z = 3.45
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = 1
coin.rotation_euler.y = 0.1
coin.rotation_euler.z = 0.1
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
'point_cache': bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print("          X                  Y                  Z                  RX                  RY                  RZ")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
else:
print("Coin is tails")


The script above does is the following:

1. Selects and deletes all extraneous objects in the scene.
2. Adds a passive cube (note that cubes are better than planes for collisions) for the ground.
3. Adds an active cylindrical coin.
• Sets friction to 0.25
• Sets restitution to 0.75
4. Adds keyframe animate to set initial velocity and rotation.
5. Prints out and tests if coin has rotated more than $$\frac{\pi}{2}$$ radians to determine heads/tails.

Driving from Mathematica

We can create a parametric model in Mathematica by replacing hard coded parameters with template variables using  delimiters as in the createCoinFlip function.

createCoinFlip[z_, rx_, ry_, rz_, friction_, restitution_] :=
StringTemplate["import bpy
from math import pi

version = bpy.app.version

for o in bpy.data.objects:
if version < (2, 80, 0):
if o.type == 'MESH' or o.type == 'EMPTY':
o.select = True
else:
o.select = False
else:
o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

if version < (2, 80, 0):
else:
bpy.ops.transform.resize(value=(1, 1, 0.1))
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"
boxObj.name = \"Ground\"

location=(0, 0, 3))
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"CYLINDER\"
bpy.context.object.rigid_body.friction = friction
bpy.context.object.rigid_body.restitution = restitution
boxObj.name = \"Coin\"
# Set reference to the coin
coin = bpy.data.objects[\"Coin\"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

sce.frame_set(4)
# Translate up a little
coin.location.z = z
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = rx
coin.rotation_euler.y = ry
coin.rotation_euler.z = rz
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
'point_cache': \
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print(\"\
X                  Y                  Z                  RX                  R\
Y                  RZ\")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
else:
print(\"Coin flip result is tails\")
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz,
"friction" -> friction, "restitution" -> restitution|>]


Blender will send a lot of information to standard out. We can parse this output with Find to extract a line of interest. Putting it all together, the following will create a python script, run Blender in the background, and parse the output.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.95, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]
(* Coin is tails *)


You can visualize the results of the simulation by removing the "--background" and repeating the step above.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.45, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --python coinflip.py >>" <> outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]


If you left click anywhere on the screen and hit the play button, then you should see the following:

Rendering in Blender

You can take advantage of Blender's photorealistic rendering capability to create a nice animation if desired.

Update: Additional Post Processing in Mathematica

Blender is geared more towards the artist whereas Mathematica is geared more towards the physicist. We can find synergy when we combine the strengths of both tools.

What follows is a simple example of how to conduct additional post-processing on a Blender simulation in Mathematica.

First, let's modify the python generation script to give the positions and orientations of the coin at each frame (we will insert a string "PosRot" to identify the appropriate lines).

createCoinFlipTransform[z_, rx_, ry_, rz_, friction_, restitution_] :=
StringTemplate["import bpy
from math import pi

version = bpy.app.version

for o in bpy.data.objects:
if version < (2, 80, 0):
if o.type == 'MESH' or o.type == 'EMPTY':
o.select = True
else:
o.select = False
else:
o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

if version < (2, 80, 0):
else:
bpy.ops.transform.resize(value=(1, 1, 0.1))
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"
boxObj.name = \"Ground\"

location=(0, 0, 3))
cylObj = bpy.context.active_object
cylObj.rigid_body.collision_shape = \"CYLINDER\"
bpy.context.object.rigid_body.friction = friction
bpy.context.object.rigid_body.restitution = restitution
cylObj.name = \"Coin\"
# Set reference to the coin
coin = bpy.data.objects[\"Coin\"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

sce.frame_set(4)
# Translate up a little
coin.location.z = z
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = rx
coin.rotation_euler.y = ry
coin.rotation_euler.z = rz
# Set Keyframes
coin.keyframe_insert(data_path=\"location\")
coin.keyframe_insert(data_path=\"rotation_euler\")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
'point_cache': \
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()

for i in range(250):
sce.frame_set(i)
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print(\"PosRot\",tr.x, tr.y, tr.z, eu.x , eu.y , eu.z )
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz,
"friction" -> friction, "restitution" -> restitution|>]


We can extract the positions and orientations of the simulation with the following code.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlipTransform[4, -Pi 0.75, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
data = ToExpression@StringSplit[#] & /@ FindList[stext, "PosRot"];
{tx, ty, tz, rx, ry, rz} = Transpose@data[[All, {2, 3, 4, 5, 6, 7}]];
Close[stext];
DeleteFile[outputfile]


We can define a cuboid and cylinder that have the same dimensions as the Blender simulation and we can create a transformation function with the following code.

box = {Cuboid[{-5, -5, -0.5}, {5, 5, 0.5}]};
cyl = {Cylinder[{{0, 0, -0.05}, {0, 0, 0.05}}, 1],
AbsolutePointSize[10],
Opacity[1], {Black, Point[{0, 0, 0}]}, {Red,
Point[{1, 0, 0}]}, {Green, Point[{0, 1, 0}]}, {Blue,
Point[{0, 0, 1}]}};
m = IdentityMatrix[4];
m[[1 ;; 3, 1 ;; 3]] = EulerMatrix[{a, b, c}, {1, 2, 3}];
m[[1 ;; 3, -1]] = {x, y, z};
transform[a_, b_, c_, x_, y_, z_] = TransformationFunction[m];


Now, we can combine plots of position and orientation (or other quantities like angular momentum) into a Manipulate[] function.

Manipulate[
Column[{Row[{ListPlot[{tx[[1 ;; i]], ty[[1 ;; i]], tz[[1 ;; i]]},
Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All,
PlotLegends -> {"tx", "ty", "tz"}],
ListPlot[{rx[[1 ;; i]], ry[[1 ;; i]], rz[[1 ;; i]]},
Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All,
PlotLegends -> {"rx", "ry", "rz"}]}],
Graphics3D[{{Opacity[0.75], Red, box},
GeometricTransformation[{Opacity[.85], Yellow, cyl},
transform[rx[[i]], ry[[i]], rz[[i]], tx[[i]], ty[[i]],
tz[[i]]]]}, SphericalRegion -> True,  Boxed -> False,
ImageSize -> {400, 400}]}], {i, 1, 250, 1}]


• Thank you very much! this answer is great! One thing: At first, running a copy of the code you posted I got an error in cmd about some "Unexpected indent". I don't know if it's something to do with how text is formatted on stackexchange, but I did fix it!
– Ivan
Apr 21, 2019 at 2:48
• I was having a little trouble getting the markdown to work here. You have SE, Mathematica, Python, and Windows all needing to escape their notion of special characters. I am glad you found it useful. I will take a look at the post tomorrow to see if I can repeat the error and try to rectify. Apr 21, 2019 at 4:00
• @Ivan I fixed the indentation problem. The "Code Sample" or ctrl-K in SE adds extra preceding spaces in the MMA code string, which causes a problem, because Python is sensitive to indentation. This is an artifact of trying to highlight a MMA code that produces a Python code in SE. Apr 21, 2019 at 13:08