Softbody physics goes haywire

118

Author

April 23, 2008 08:43 AM

Hi guys, I found this funny little demo app that shows how to implement 2D soft body physics: link. Seeing as how I'm trying to get comfortable with Python and Pyglet, I tried to port the program to Python. Thing is, the ball now blows up after a little while :( I've tried to compare my implementation with the original, but I just can't find what's going wrong. I could really use some help with this.


from pyglet.gl import *
import pyglet

from math import pi, sin, cos, sqrt, pow, fabs

#config = Config(double_buffer=True,)


# constants
screenSizeFac = 7
numPoints     = 20
mass          = 1.0
ballRadius    = 0.5
kS            = 755.0
kD            = 35.0
gy            = -10.0
dt            = 0.005
finalPressure = 45.0

pressure = 0.0

winWidth  = 500
winHeight = 500


window = pyglet.window.Window(winWidth,winHeight, resizable=True, vsync=False)

clock = pyglet.clock.ClockDisplay()


class Point2D:
  def __init__(self, x=0.0, y=0.0, vx=0.0, vy=0.0, fx=0.0, fy=0.0):
    self.x = x
    self.y = y
    self.oldX = x
    self.oldY = y
    self.vx = vx
    self.vy = vy
    self.fx = fx
    self.fy = fy
    
  def __str__(self):
    return '(%d,%d)' % (self.x, self.y)
    
class Spring:
  def __init__(self, i=0, j=0, length=0.0, nx=0.0, ny=0.0):
    self.i = i
    self.j = j
    self.length = length
    self.nx = nx
    self.ny = ny

  def __str__(self):
    return '(%d,%d,%f)' % (self.i, self.j, self.length)

mySprings = []
myPoints  = []
centerPoint = Point2D()

def addSpring(i, j):
  # compute spring length from distance between points
  length = sqrt( (myPoints.x - myPoints[j].x)*(myPoints.x - myPoints[j].x) + 
                 (myPoints.y - myPoints[j].y)*(myPoints.y - myPoints[j].y) )

  # connect the points with springs
  mySprings.append( Spring(i, j, length) )


def createBall():
  # generate the soft body
  for i in range(0, numPoints):
    
    # generate point
    pointX = ballRadius*sin( i*(2*pi)/(numPoints-1) )
    pointY = ballRadius*cos( i*(2*pi)/(numPoints-1) ) + screenSizeFac / 2
    myPoints.append( Point2D(pointX, pointY) )
  
  for i in range(0, numPoints-1):
    addSpring(i,i+1)

  #final spring (between last and first point)
  addSpring(numPoints-1, 1)
  
  for i in range(0, numPoints):
    print mySprings
  
  print myPoints[2]
  
  #create center point for rendering
  centerPoint.x = myPoints[0].x + (myPoints[numPoints/2-1].x-myPoints[0].x)*0.5
  centerPoint.y = myPoints[0].y + (myPoints[numPoints/2-1].y-myPoints[0].y)*0.5
  
## Calculate all the forces acting on the ball
def accumulateForces():
  r12d = 0.0 #length of p1 - p2 vector
  volume = 0.0      #volume of the body
  pressurev = 0.0   #pressure force value

  # gravity  
  for i in range(0, numPoints):
    if (pressure - finalPressure &gt;= 0):
      myPoints.fy = mass*gy
    else:
      myPoints.fy = 0

  # spring force
  for spring in mySprings:
    x1 = myPoints[spring.i].x
    y1 = myPoints[spring.i].y
    x2 = myPoints[spring.j].x
    y2 = myPoints[spring.j].y
    
    r12d = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) )
    
    if (r12d != 0):
      vx12 = myPoints[spring.i].vx - myPoints[spring.j].vx
      vy12 = myPoints[spring.i].vy - myPoints[spring.j].vy
      
      f = (r12d - spring.length)*kS + (vx12*(x1-x2) + vy12*(y1-y2)) * kD/r12d

      Fx = ((x1-x2) / r12d) * f
      Fy = ((y1-y2) / r12d) * f
      
      myPoints[spring.i].fx -= Fx
      myPoints[spring.i].fy -= Fy
      myPoints[spring.j].fx += Fx
      myPoints[spring.j].fy += Fy
    
    spring.nx = (y1-y2) / r12d
    spring.ny = -(x1-x2) / r12d
  
  #pressure force
  for spring in mySprings:
    # calculate total volume
    x1 = myPoints[spring.i].x
    y1 = myPoints[spring.i].y
    x2 = myPoints[spring.j].x
    y2 = myPoints[spring.j].y
 
    r12d = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) )
    
    volume += 0.5 * fabs(x1-x2) * fabs(spring.nx) * r12d
    
  for spring in mySprings:
    x1 = myPoints[spring.i].x
    y1 = myPoints[spring.i].y
    x2 = myPoints[spring.j].x
    y2 = myPoints[spring.j].y
 
    r12d = sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) )

    pressurev = r12d * pressure * (1.0/volume)
    myPoints[spring.i].fx += spring.nx * pressurev
    myPoints[spring.i].fy += spring.ny * pressurev
    myPoints[spring.j].fx += spring.nx * pressurev
    myPoints[spring.j].fy += spring.ny * pressurev
  
 
 
## Do integration
def euler():
  drx = 0.0
  dry = 0.0

  for point in myPoints:
    # calc x
    point.vx += (point.fx / mass)*dt
    drx = point.vx*dt

    #drx = 2*point.x - point.oldX + (point.fx / mass)
    #point.oldX = point.x
    
    # check border x
    if (point.x + drx &lt; -screenSizeFac):
      drx = -screenSizeFac - point.x
      point.vx = -0.1*point.vx
      point.vy = 0.95*point.vy
    elif (point.x + drx &gt; screenSizeFac):
      drx = screenSizeFac - point.x
      point.vx = -0.1*point.vx
      point.vy = 0.95*point.vy

    point.x += drx

    # do the same for y
    point.vy += (point.fy / mass)*dt
    dry = point.vy*dt

    #dry = 2*point.y - point.oldY + (point.fy / mass)
    #point.oldY = point.y
    
    # check border y
    if (point.y + drx &lt; -screenSizeFac):
      dry = -screenSizeFac - point.y
      point.vx = 0.95*point.vx
      point.vy = -0.1*point.vy
    elif (point.x + drx &gt; screenSizeFac):
      dry = screenSizeFac - point.y
      point.vx = 0.95*point.vx
      point.vy = -0.1*point.vy

    point.y += dry
    
    if (point.x &gt; screenSizeFac):
      point.x = screenSizeFac
    if (point.x &lt; -screenSizeFac):
      point.x = -screenSizeFac
    if (point.y &gt; screenSizeFac):
      point.y = screenSizeFac
    if (point.y &lt; -screenSizeFac):
      point.y = -screenSizeFac

def updatePhysics(dt):
  global pressure, finalPressure, centerPoint
  accumulateForces()
  euler()

  # update middle point
  centerPoint.x = myPoints[0].x + (myPoints[numPoints/2-1].x-myPoints[0].x)*0.5
  centerPoint.y = myPoints[0].y + (myPoints[numPoints/2-1].y-myPoints[0].y)*0.5

  if(pressure &lt; finalPressure):
    pressure += finalPressure/300.0
    print 'pressure is now: ',pressure

@window.event
def on_resize(width, height):
  # Override the default on_resize handler to create a 3D projection
  glViewport(0, 0, width, height)
  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()
  gluOrtho2D(-screenSizeFac, screenSizeFac, -screenSizeFac, screenSizeFac)
  glMatrixMode(GL_MODELVIEW)

  return pyglet.event.EVENT_HANDLED


@window.event
def on_draw():
  glClearColor(0,0,0,0)
  glClear(GL_COLOR_BUFFER_BIT)

  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()
  gluOrtho2D(0, winWidth, 0, winHeight)
  glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)

  clock.draw()

  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()
  gluOrtho2D(-screenSizeFac, screenSizeFac, -screenSizeFac, screenSizeFac)

  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
  glBegin(GL_TRIANGLES)
  
  for i in range(0, numPoints):
    glColor3f(0.8,0.4,0.4)
    
    glVertex2f( myPoints[mySprings.i].x, myPoints[mySprings.i].y )
    glVertex2f( myPoints[mySprings.j].x, myPoints[mySprings.j].y )
    glVertex2f( centerPoint.x, centerPoint.y )
    
    #glVertex2f( myPoints[numPoints-mySprings.i+1].x, myPoints[numPoints-mySprings.i+1].y )
    #glVertex2f( myPoints[numPoints-mySprings.j+1].x, myPoints[numPoints-mySprings.j+1].y )
  
  glEnd()
  
createBall()

pyglet.clock.schedule_interval(updatePhysics, 1.0/60)

pyglet.app.run()

Antheus

2,410

April 23, 2008 01:45 PM

It's just a case of instability during integration. There are some ways to fix this reliably, but rarely for complex systems. Realistically, you have several options:

- Change the integration method.
- If using variable time-step, use fixed time step
- Improve precision of your floating point numbers
- Use fixed-point math
- Reduce the time step
- Add more dampening
- Put artificial limits on parameters of your system
- Monitor your system for rapidly growing values

With numerical integration, systems explode. As said, mathematically it is sometimes possible to determine maximum safe time step. But after taking into consideration the inaccuracy of floating point numbers, as well as general complexity of such systems, it's often not viable.

For almost exactly the same type of simulation, I used double precision numbers, used RK4, added some artificial dampening, and, since it was possible in my case due to domain-specific knowledge, I monitored the change in acceleration values. If those started growing more than a certain limit, I'd cut the time-step down drastically. If they were stable, I'd increase the time step.

Menno

118

Author

April 23, 2008 05:16 PM

Thanks for the explanation. I just thought it was weird that my version exploded, whilst the example program I based it on runs fine. Perhaps its a precision difference between Python and C?

Antheus

2,410

April 23, 2008 08:18 PM

Possible, or perhaps you have a bug.

See here for an explanation of instability itself. It really isn't an easy problem.

Try reducing the step to one half of what you have now. Perhaps try multiplying the velocities on each iteration by .99 to add a bit more dampening if you don't already do that.

There's usually some tweaking involved with such systems.

Softbody physics goes haywire

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Softbody physics goes haywire

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines