Skip to main content

packages - 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!



Answer




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.



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

for o in bpy.data.objects:
    if o.type == 'MESH' or o.type == 'EMPTY':

        o.select = True
    else:
        o.select = False

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

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))

bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
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')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)
# 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.:
    print("Coin is heads")
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.



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

for o in bpy.data.objects:

    if o.type == 'MESH' or o.type == 'EMPTY':
        o.select = True
    else:
        o.select = False

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

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))

bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"
boxObj.name = \"Ground\"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, \
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')

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')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)

# 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.:
    print(\"Coin flip result is heads\")
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];
stext = OpenRead[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];
stext = OpenRead[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:


Blender Coin Flip



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


Rendered Coin Flip



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

for o in bpy.data.objects:
    if o.type == 'MESH' or o.type == 'EMPTY':
        o.select = True
    else:
        o.select = False


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

# Add the floor
bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = \"BOX\"

boxObj.name = \"Ground\"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, \
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
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')

# Advance two frames and add translational and rotational motion
sce.frame_set(3)
# 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];
stext = OpenRead[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}]

Physics Graphs and 3D Animation



Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...
Post a Comment
Read more