ClothObject.py
from OpenGL.GLUT import *
from OpenGL.GL import *
from OpenGL.GLU import *
import numpy as np
import math
################### new functions for Lab 4 ################
def convertLongVector(x) :
# concatenate all 3-d vectors stored in array x
# to create a large vector longX
longX = np.zeros(shape=(len(x)*3,))
for i in range(0, len(x)) :
longX[i*3:i*3+3] = x[i]
return longX
def set33SubMatrix(M, i, j, m) :
# set the submatrix of a larg matrix M at the i-th row and the j-th column
# with a 3x3 small matrix m
for row in range(0,3) :
for col in range(0,3) :
M[i*3+row, j*3+col] = m[row,col]
def set33SubMatrixSymmetric(M, i, j, m):
set33SubMatrix(M, i, j, m)
set33SubMatrix(M, j, i, m)
#############################################################
class ClothObject :
def __init__(self, sizeX, sizeZ, nW, nH):
self.dx = sizeX / (nW-1)
self.dz = sizeZ / (nH-1)
self.dDiag = math.sqrt(self.dx**2.0 + self.dz**2.0)
self.damp = 0.0
self.nW, self.nH = nW, nH
self.clothSize = (sizeX, sizeZ)
# vertex array (x,y,z)
self.verts = np.zeros(shape=(nW * nH, 3))
self.velocity = np.zeros(shape=(nW*nH, 3))
self.force = np.zeros(shape=(nW*nH, 3))
self.totalMass = 1.0
self.stiffness = 1.0
self.totalNumberOfParticles = nW*nH
self.mass = np.zeros(shape=(len(self.verts),))
particleMass = self.totalMass / self.totalNumberOfParticles
for i in range(0, self.totalNumberOfParticles) :
self.mass[i] = particleMass
# spring links (i, j)
self.springs = np.zeros(shape=((nW - 1) * nH + (nH - 1) * nW + (nW - 1) * (nH - 1) * 2, 2), dtype=np.int32)
# relaxed length of spring : l0
self.l0 = np.zeros(shape=(len(self.springs),))
self.numXSprings = (self.nW - 1) * self.nH
self.numZSprings = (self.nH - 1) * self.nW
self.numDiagonalSprings = (self.nW - 1) * (self.nH - 1)
self.resetMassSpring()
def resetMassSpring(self):
for i in range(0, self.nW * self.nH):
self.verts[i] = [(i % self.nW) * self.dx, 1, (int(i / self.nW)) * self.dz]
self.velocity[i] = [0,0,0]
for i in range(0, self.numXSprings):
row = int(i / (self.nW - 1))
col = i % (self.nW - 1)
self.springs[i] = [row * self.nW + col, row * self.nW + col + 1]
self.l0[i] = self.dx
for i in range(0, self.numZSprings):
row = int(i / (self.nW))
col = i % (self.nW)
self.springs[self.numXSprings + i] = [row * self.nW + col, (row + 1) * self.nW + col]
self.l0[(self.nW - 1) * self.nH + i] = self.dz
for i in range(0, self.numDiagonalSprings):
row = int(i / (self.nW - 1))
col = i % (self.nW - 1)
self.springs[self.numXSprings + self.numZSprings + 2 * i] = [row * self.nW + col,
(row + 1) * self.nW + col + 1]
self.l0[self.numXSprings + self.numZSprings + 2 * i] = self.dDiag
self.springs[self.numXSprings + self.numZSprings + 2 * i + 1] = [row * self.nW + col + 1,
(row + 1) * self.nW + col]
self.l0[self.numXSprings + self.numZSprings + 2 * i + 1] = self.dDiag
def drawSpring(self):
glColor3f(1, 1, 1)
glBegin(GL_LINES)
for i in range(0, len(self.springs)) :
idx0 = self.springs[i][0]
idx1 = self.springs[i][1]
loc0 = self.verts[idx0]
loc1 = self.verts[idx1]
glVertex3fv(loc0)
glVertex3fv(loc1)
glEnd()
def computeForce(self):
## reset forces
for i in range(0, self.totalNumberOfParticles) :
self.force[i] = np.array([0.0, 0.0, 0.0])
## compute spring forces
for i in range(0, len(self.springs)) :
idx0 = self.springs[i][0]
idx1 = self.springs[i][1]
xi = self.verts[idx0]
xj = self.verts[idx1]
xji = xj - xi
vi = self.velocity[idx0]
vj = self.velocity[idx1]
vji = vj - vi
l = np.linalg.norm(xji)
dir = xji / l
deform = l - self.l0[i]
f = self.stiffness * deform * dir + self.damp * vji
self.force[idx0] = self.force[idx0] + f
self.force[idx1] = self.force[idx1] - f
################### new function for Lab 6 ################
def computePhi(self,dt):
phi = np.zeros(shape=(self.totalNumberOfParticles,3), dtype=np.float64)
for e in range(0, len(self.springs)):
i, j = self.springs[e][0], self.springs[e][1]
xji = self.verts[j] - self.verts[i]
vji = self.velocity[j] - self.velocity[i]
lt = np.linalg.norm(xji)
mi, mj = self.mass[i], self.mass[j]
M = (mi*mj) / (mi+mj)
d = lt - self.l0[e]
xjiHat = xji / lt
dldt = np.dot(xjiHat, vji)
alpha = self.damp / (2.0*M)
omega = math.sqrt(self.stiffness / M - alpha**2.0)
A = math.sqrt(d*d + (dldt*dldt)*M / self.stiffness)
if A != 0.0: dDA = d / A
else: dDA = 0.0
arcsine = math.asin(dDA)
if dldt < 0.0: arcsine = np.pi - arcsine
elif d < 0.0 and dldt >= 0.0: arcsine = 2.0 * np.pi + arcsine
T0 = arcsine / omega
coeff = -A * self.stiffness / (omega ** 2.0 + alpha ** 2.0)
Et, Eh = math.exp(-alpha * T0), math.exp(-alpha * dt)
st, ct = math.sin(omega * T0), math.cos(omega * T0)
sh, ch = math.sin(omega * dt), math.cos(omega * dt)
phi_ij = coeff * Et * (
st * (Eh * (alpha * ch - omega * sh) - alpha) +
ct * (Eh * (alpha * sh + omega * ch) - omega));
phi[i] += (phi_ij * mj/(mi+mj)) * xjiHat
phi[j] -= (phi_ij * mi/(mi+mj)) * xjiHat
return phi
def integrateHarmonicOscillation(self, dt) :
dVel = np.zeros(shape=(self.totalNumberOfParticles * 3), dtype=np.float64)
phi = self.computePhi(dt)
for i in range(self.totalNumberOfParticles):
dVel[i * 3:i * 3 + 3] = phi[i] / self.mass[i]
return dVel
#############################################################
def update(self, dt):
gravity = np.array([0.0, -9.8, 0.0])
self.computeForce()
################### new function call for Lab 6 ################
dVel = self.integrateHarmonicOscillation(dt)
#############################################################
gravity = np.array([0.0, -9.8, 0.0])
for i in range (0, self.totalNumberOfParticles) :
if i is not 0 and i is not self.nW-1 : # these two masses are static
self.velocity[i] = self.velocity[i] + dVel[i*3:i*3+3] + gravity*dt
self.verts[i] = self.verts[i] + self.velocity[i]*dt
#### Simple Euler method for gravity simulation
main.py
import wx # requires wxPython package
from wx import glcanvas
from OpenGL.GL import *
from OpenGL.GLU import *
import ClothObject
class MyFrame(wx.Frame) :
def __init__(self):
self.size = (1280, 720)
wx.Frame.__init__(self, None, title = "wx frame", size = self.size,
style = wx.DEFAULT_FRAME_STYLE ^ wx.RESIZE_BORDER)
self.panel = MyPanel(self)
class MyPanel(wx.Panel) :
def __init__(self, parent):
wx.Panel.__init__(self, parent)
self.canvas = OpenGLCanvas(self)
self.bAnimation = False
self.resetButton = wx.Button(self, wx.ID_ANY, "Mass-Spring Reset", pos=(1030, 20), size=(200,40), style=0 )
self.animationButton = wx.Button(self, wx.ID_ANY, "Animate/Stop", pos=(1030, 60), size=(200, 40), style=0)
self.stiffnessLabel = wx.StaticText(self, -1, pos=(1030, 150), style=wx.ALIGN_CENTER)
self.stiffnessLabel.SetLabel("Stiffness: " + str(self.canvas.clothObject.stiffness))
self.stiffnessSlider = wx.Slider(self, -1, pos=(1030, 180), size=(200,50),
style = wx.SL_HORIZONTAL|wx.SL_AUTOTICKS,
value=1, minValue=1, maxValue = 30)
self.stepLabel = wx.StaticText(self, -1, pos=(1030, 250), style=wx.ALIGN_CENTER)
self.stepLabel.SetLabel("Time Interval: " + str(self.canvas.stepSize))
self.stepSlider = wx.Slider(self, -1, pos=(1030, 280), size=(200, 50),
style=wx.SL_HORIZONTAL | wx.SL_AUTOTICKS,
value=10, minValue=1, maxValue=300)
self.dampLabel = wx.StaticText(self, -1, pos=(1030, 350), style=wx.ALIGN_CENTER)
self.dampLabel.SetLabel("Damping: " + str(self.canvas.clothObject.damp))
self.dampSlider = wx.Slider(self, -1, pos=(1030, 380), size=(200, 50),
style=wx.SL_HORIZONTAL | wx.SL_AUTOTICKS,
value=0, minValue=0, maxValue=30)
self.Bind(wx.EVT_BUTTON, self.OnAnimationButton, self.animationButton)
self.Bind(wx.EVT_BUTTON, self.OnResetButton, self.resetButton)
self.Bind(wx.EVT_SLIDER, self.OnStiffnessSlider, self.stiffnessSlider)
self.Bind(wx.EVT_SLIDER, self.OnStepSlider, self.stepSlider)
self.Bind(wx.EVT_SLIDER, self.OnDampSlider, self.dampSlider)
def OnAnimationButton(self, event) :
if self.bAnimation is False :
self.bAnimation = True
else :
self.bAnimation = False
self.canvas.bAnimation = self.bAnimation
def OnResetButton(self, event) :
self.canvas.clothObject.resetMassSpring()
def OnStiffnessSlider(self, event):
val = event.GetEventObject().GetValue()
stiffness = 2**val
self.stiffnessLabel.SetLabel("Stiffness :" + str(stiffness))
self.canvas.clothObject.stiffness = stiffness
def OnStepSlider(self, event):
val = event.GetEventObject().GetValue()
stepSize = 0.0001*val
self.stepLabel.SetLabel("Time Interval :" + str(stepSize))
self.canvas.stepSize = stepSize
def OnDampSlider(self, event):
val = event.GetEventObject().GetValue()
damp = val/10.0
self.dampLabel.SetLabel("Damping Coeff. :" + str(damp))
self.canvas.clothObject.damp = damp
class OpenGLCanvas(glcanvas.GLCanvas):
def __init__(self, parent) :
self.initialized = False
self.size = (1024,720)
self.aspect_ratio = 1
self.stepSize = 0.001
glcanvas.GLCanvas.__init__(self, parent, -1, size = self.size)
self.context = glcanvas.GLContext(self)
self.SetCurrent(self.context)
self.Bind(wx.EVT_PAINT, self.OnDraw)
self.Bind(wx.EVT_IDLE, self.OnIdle)
self.InitGL()
self.clothObject = ClothObject.ClothObject(1.0,1.0,10, 10)
self.bAnimation = False
def setAnimationFlag(self, bOption):
self.bAnimation = bOption
def resetMesh(self):
self.clothObject.resetMassSpring()
def InitGL(self):
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
self.aspect_ratio = float(self.size[0]) / self.size[1]
gluPerspective(60, self.aspect_ratio, 0.1, 100.0)
glViewport(0,0,self.size[0], self.size[1])
def OnDraw(self, event):
# clear color and depth buffers
if not self.initialized :
self.InitGL()
self.initialized = True
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
# position viewer
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
gluLookAt(2,2,2, 0,0,0, 0,1,0)
self.clothObject.drawSpring()
self.SwapBuffers()
def OnIdle(self, event):
if self.bAnimation :
self.clothObject.update(self.stepSize)
self.Refresh()
def main() :
app = wx.App()
frame = MyFrame()
frame.Show()
app.MainLoop()
if __name__ == "__main__" :
main()