意味も分からず線とか追加。厳密な当たり判定すらできず、ズルしてみたけど、うまく動いとる。ゲームとかならまだしも、解析でやっちゃいかんでしょー。
# -*- coding: utf-8 -*-
import sys
import math
import random
from PySide import QtGui,QtCore
g = 9.80665 # 重力加速度
rho = 10 # 粒子間密度
pn = 0 #粒子数
class Particle():
def __init__(self,x,y,vx=0,vy=0):
global pn
self.n = pn #粒子NO.
pn += 1
self.r = 5 #半径
self.x = x #X座標
self.y = y #Y座標
self.a = 0 #角度
self.dx = 0 #X方向増加量
self.dy = 0 #Y方向増加量
self.da = 0 #角度増加量
self.vx = vx #X方向速度
self.vy = vy #Y方向速度
self.va = 0 #角速度
self.fy = 0
self.fx = 0
self.fm = 0
self.m = 4.0/3.0*math.pi*rho*self.r**3 # 質量
self.Ir = math.pi*rho*self.r**4/2.0 #慣性モーメント
self.en = [] #弾性力（直方向）
self.es = [] #弾性力（せん断方向）
def config(self):
self.en = [0 for i in range(pn)]
self.es = [0 for i in range(pn)]
def nextStep(self,dt):
#位置更新（オイラー差分）
ax = self.fx/self.m
ay = self.fy/self.m
aa = self.fm/self.Ir
self.vx += ax*dt
self.vy += ay*dt
self.va += aa*dt
self.dx = self.vx*dt
self.dy = self.vy*dt
self.da = self.va*dt
self.x += self.dx
self.y += self.dy
self.a += self.da
class Line():
def __init__(self,x1,y1,x2,y2):
global pn
self.n = pn #粒子NO.
pn += 1
self.x1 = x1
self.y1 = y1
self.x2 = x2
self.y2 = y2
class Interface():
def __init__(self):
pass
class DEM():
def __init__(self):
self.step = 0
self.dt = 0.001
self.pars = []
self.lines = []
# set particle
for i in range(10):
p = Particle(random.randint(-140,140),random.randint(0,90),
random.randint(-40,40),random.randint(-40,40))
self.pars.append(p)
p1 = Particle(120,80,0,0)
p2 = Particle(-40,50,0,0)
#p3 = Particle(10,20,-20,30)
#p4 = Particle(-80,80,40,5)
#self.pars.append(p1)
self.pars.append(p2)
#self.pars.append(p3)
#self.pars.append(p4)
# set line
l1 = Line(-150,-100, 150,-100)
l2 = Line(-150, 100,-150,-100)
l3 = Line( 150, 100, 150,-100)
l4 = Line(-150, 100, 150, 100)
l5 = Line(-120,-10,80,-40)
l6 = Line(100,50,-30,30)
self.lines.append(l1)
self.lines.append(l2)
self.lines.append(l3)
self.lines.append(l4)
self.lines.append(l5)
self.lines.append(l6)
for p in self.pars:
p.config()
def calcForce(self):
#2粒子間の接触判定
combs = [(p1, p2) for p1 in self.pars for p2 in self.pars]
for p1,p2 in combs:
if p1.n == p2.n:
continue
lx = p1.x - p2.x
ly = p1.y - p2.y
ld = math.sqrt(lx**2+ly**2)
if (p1.r+p2.r)>ld:
cos_a = lx/ld
sin_a = ly/ld
self.force2par(p1,p2,cos_a,sin_a)
else:
p1.en[p2.n] = 0.0
p1.es[p2.n] = 0.0
#粒子と線の接触判定
combs = [(l, p) for l in self.lines for p in self.pars]
for l,p in combs:
hit = False
th0 = math.atan2(l.y2-l.y1,l.x2-l.x1)
th1 = math.atan2(p.y-l.y1,p.x-l.x1)
a = math.sqrt((p.x-l.x1)**2+(p.y-l.y1)**2)
d = abs(a*math.sin(th1-th0))
if d < p.r:
b = math.sqrt((p.x-l.x2)**2+(p.y-l.y2)**2)
s = math.sqrt((l.x2-l.x1)**2+(l.y2-l.y1)**2)
if a < s and b < s:
s1 = math.sqrt(a**2-d**2)
x = l.x1 + s1*math.cos(th0)
y = l.y1 + s1*math.sin(th0)
hit = True
elif a < b and a < p.r:
x = l.x1
y = l.y1
hit = True
elif b < p.r:
x = l.x2
y = l.y2
hit = True
if hit:
self.force2line(p,l,x,y)
else:
p.en[l.n] = 0.0
p.es[l.n] = 0.0
pass
#外力
for p in self.pars:
p.fy += -g*p.m #重力
def force2line(self,p1,line,x,y):
kn = 10**7 #弾性係数（法線方向）
etan = 0 #粘性係数（法線方向）
ks = 50000 #弾性係数（せん断方向）
etas = 1000 #粘性係数（せん断方向）
frc = 10 #摩擦係数
lx = p1.x - x
ly = p1.y - y
ld = math.sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
#相対的変位増分
un = +(p1.dx)*cos_a+(p1.dy)*sin_a
us = -(p1.dx)*sin_a+(p1.dy)*cos_a+(p1.r*p1.da)
#相対的速度増分
vn = +(p1.vx)*cos_a+(p1.vy)*sin_a
vs = -(p1.vx)*sin_a+(p1.vy)*cos_a+(p1.r*p1.va)
#合力
p1.en[line.n] += kn*un
p1.es[line.n] += ks*us
hn = p1.en[line.n] + etan*vn
hs = p1.es[line.n] + etas*vs
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif abs(hs)-frc*hn >= 0.0:
#摩擦力以上のせん断力は働かない
hs = frc*abs(hn)*hs/abs(hs)
#合力
p1.fx += -hn*cos_a + hs*sin_a
p1.fy += -hn*sin_a - hs*cos_a
p1.fm -= p1.r*hs
def force2par(self,p1,p2,cos_a,sin_a):
kn = 10**6 #弾性係数（法線方向）
etan = 1000 #粘性係数（法線方向）
ks = 5000 #弾性係数（せん断方向）
etas = 1000 #粘性係数（せん断方向）
frc = 10 #摩擦係数
#相対的変位増分
un = +(p1.dx-p2.dx)*cos_a+(p1.dy-p2.dy)*sin_a
us = -(p1.dx-p2.dx)*sin_a+(p1.dy-p2.dy)*cos_a+(p1.r*p1.da+p2.r*p2.da)
#相対的速度増分
vn = +(p1.vx-p2.vx)*cos_a+(p1.vy-p2.vy)*sin_a
vs = -(p1.vx-p2.vx)*sin_a+(p1.vy-p2.vy)*cos_a+(p1.r*p1.va+p2.r*p2.va)
#合力
p1.en[p2.n] += kn*un
p1.es[p2.n] += ks*us
hn = p1.en[p2.n] + etan*vn
hs = p1.es[p2.n] + etas*vs
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif abs(hs)-frc*hn >= 0.0:
#摩擦力以上のせん断力は働かない
hs = frc*abs(hn)*hs/abs(hs)
#合力
p1.fx += -hn*cos_a + hs*sin_a
p1.fy += -hn*sin_a - hs*cos_a
p1.fm -= p1.r*hs
def clacStep(self):
for p in self.pars:
p.fx = 0
p.fy = 0
p.fm = 0
self.calcForce()
for p in self.pars:
p.nextStep(self.dt)
self.step += 1
class MainWindow(QtGui.QMainWindow):
def __init__(self, *argv, **keywords ):
super(MainWindow,self).__init__(*argv,**keywords)
self.resize(320, 220)
self.view = QtGui.QGraphicsView(self)
self.scene = QtGui.QGraphicsScene(self.view)
self.view.setScene(self.scene)
self.setCentralWidget(self.view)
self.pen = QtGui.QPen()
self.pen.setColor(QtGui.QColor('#000000'))
self.dem = DEM()
self.step_limit = int(1/self.dem.dt*5)
#self.step_limit = float("inf")
print('limit',self.step_limit)
self.output_step = 0
self.output_interval = 10
self.redraw()
self.timer = QtCore.QTimer(self)
self.timer.setInterval(int(1000*self.dem.dt))
print('interval',int(1000*self.dem.dt))
self.timer.timeout.connect(self.drawLoop)
self.timer.start()
def redraw(self):
for item in self.scene.items():
self.scene.removeItem(item)
pass
for p in self.dem.pars:
self.scene.addEllipse(p.x-p.r,-p.y-p.r,p.r*2,p.r*2,self.pen)
for l in self.dem.lines:
self.scene.addLine(l.x1,-l.y1,l.x2,-l.y2,self.pen)
def drawLoop(self):
self.dem.clacStep()
self.redraw()
if self.dem.step % self.output_interval == 0:
#self.saveJpeg()
self.output_step += 1
if self.dem.step >= self.step_limit:
self.timer.stop()
self.close()
def saveJpeg(self):
pix = QtGui.QPixmap(330,220)
pix.fill(QtCore.Qt.white)
paint = QtGui.QPainter(pix)
self.scene.render(paint)#,QtCore.QRect(0,0,300,200))
pix.save('c://Temp/dem/capture%03d.jpg' % self.output_step,'JPG')
paint.end()
def main():
app = QtGui.QApplication(sys.argv)
w = MainWindow()
w.show()
sys.exit(app.exec_())
if __name__ == '__main__':
main()