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

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

equation solving - Invert and fit implicitly defined curve

I need to fit an implicitly defined curve. I thought I could get some data out of Solve , and then using FindFit . Therefore, I would like to find the relation the parametric curve defined by $F(x,y)=0$: Solve[-(1/2) + 1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + (0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0, y] But I can't get an output: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >> Edit: In particular, I would like to fit the data coming from the curve with the expression of another curve, and not with a function $f(x)$. In particular, since this clearly looks like a cardioid , I would like it to fit to something like it. What other strategies could I try?
Post a Comment
Read more

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...
Post a Comment
Read more