３といっても新しいことはしていないのだけど。
前回記事のPythonプログラムは、PySideというGUIモジュールを使っているのだけど、考えてみれば、解析の様子を見せるだけなので、標準GUIの方が良いよなぁっと思って、Tkinterに書き換えました。・・・それだけ、です。
PureなPython(2.7)で動くプログラムになりました。
# -*- coding: utf-8 -*-
import sys
import math
import random
import time
import Tkinter
G = 9.80665 # 重力加速度
rho = 10 # 粒子間密度
elem_num = 0 #粒子数
class Element(object):
def __init__(self):
global elem_num
self.n = elem_num #要素No.
elem_num += 1
self.r = 0 #半径
self.x = 0 #X座標
self.y = 0 #Y座標
self.a = 0 #角度
self.dx = 0 #X方向増加量
self.dy = 0 #Y方向増加量
self.da = 0 #角度増加量
self.vx = 0 #X方向速度
self.vy = 0 #Y方向速度
self.va = 0 #角速度
self.fy = 0
self.fx = 0
self.fm = 0
self.en = [] #弾性力（直方向）
self.es = [] #弾性力（せん断方向）
class Particle(Element):
def __init__(self,x,y,vx=0,vy=0):
super(Particle,self).__init__()
self.type = 1
self.x = x #X座標
self.y = y #Y座標
self.vx = vx #X方向速度
self.vy = vy #Y方向速度
self.r = 5 #半径
self.m = 4.0/3.0*math.pi*rho*self.r**3 # 質量
self.Ir = math.pi*rho*self.r**4/2.0 #慣性モーメント
def config(self):
self.en = [0 for i in range(elem_num)]
self.es = [0 for i in range(elem_num)]
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(Element):
def __init__(self,x1,y1,x2,y2):
super(Line,self).__init__()
self.type = 2
self.x1 = x1
self.y1 = y1
self.x2 = x2
self.y2 = y2
class Interface():
def set(self,elem1,elem2):
if elem1.type == 1 and elem2.type == 1:
#粒子間
self.kn = 10**6 #弾性係数（法線方向）
self.etan = 10000 #粘性係数（法線方向）
self.ks = 10*6 #弾性係数（せん断方向）
self.etas = 10000 #粘性係数（せん断方向）
self.frc = 10*6 #摩擦係数
elif elem1.type == 1 and elem2.type == 2:
#粒子-Line間
self.kn = 10**7
self.etan = 0
self.ks = 10*6
self.etas = 0
self.frc = 10*6
class DEM():
def __init__(self):
self.step = 0
self.dt = 0.001
self.pars = []
self.lines = []
self.beams = []
self.interface = Interface()
# set particle random
for i in range(10):
p = Particle(random.randint(20,180),random.randint(20,180),
random.randint(-40,40),random.randint(-40,40))
self.pars.append(p)
# set line
l1 = Line(5,5, 295,5)
l2 = Line(5,195,5,5)
l3 = Line(295,195,295,5)
l4 = Line(5,195,295,195)
l5 = Line(10,150,180,140)
l6 = Line(30,30,160,50)
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()
return
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:
lx = p.x - x
ly = p.y - y
ld = math.sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
self.force2par(p,l,cos_a,sin_a)
else:
p.en[l.n] = 0.0
p.es[l.n] = 0.0
#外力
for p in self.pars:
p.fy += -G*p.m #重力
for b in self.beams:
b.fy += -G*b.m
def force2par(self,p1,p2,cos_a,sin_a):
#相対的変位増分
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)
self.interface.set(p1,p2)
#合力（局所座標系）
p1.en[p2.n] += self.interface.kn*un
p1.es[p2.n] += self.interface.ks*us
hn = p1.en[p2.n] + self.interface.etan*vn
hs = p1.es[p2.n] + self.interface.etas*vs
if hn <= 0.0: #法線力がなければ、せん断力は０ hs = 0.0 elif abs(hs) >= self.interface.frc*hn:
#摩擦力以上のせん断力は働かない
hs = self.interface.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(Tkinter.Tk):
def __init__(self):
Tkinter.Tk.__init__(self)
self.dem = DEM()
self.canvas = Tkinter.Canvas(self, bg="white")
self.canvas.pack(fill=Tkinter.BOTH,expand=True)
self.geometry('300x200')
h = 200
for l in self.dem.lines:
self.canvas.create_line(l.x1,h-l.y1,l.x2,h-l.y2,width=1)
self.redraw()
self.st = time.time()
print('init')
def redraw(self):
self.canvas.delete('elem')
h = 200
for p in self.dem.pars:
x1 = round(p.x-p.r)
y1 = round(h-p.y-p.r)
x2 = round(p.x+p.r)
y2 = round(h-p.y+p.r)
self.canvas.create_oval(x1,y1,x2,y2,tags='elem')
x1 = round(p.x)
y1 = round(h-p.y)
x2 = round(p.x+p.r*math.cos(p.a))
y2 = round(h-p.y+p.r*math.sin(p.a))
self.canvas.create_line(x1,y1,x2,y2,tags='elem')
self.update_idletasks()
self.after(100,self.redraw)
def calcloop(self):
self.dem.clacStep()
if self.dem.step % 1000 == 0:
print(self.dem.step)
#print('time:{0}'.format(time.time()-self.st))
self.st = time.time()
self.after(0,self.calcloop)
def main():
w = MainWindow()
w.after(500,w.calcloop)
w.after(500,w.redraw)
w.mainloop()
if __name__ == '__main__':
main()
※円形要素はランダムな位置に作成されるので、たまに位置がカブってエラーになる時がある
Tkinterにすると一つ困るのが、画像の出力。描いたものを直接PNG等に出来ないので、お手軽に動画にできない。うーむ。
追記：お手軽な解決方法がありました。スクリーンショットを撮るっていう。ピュアPythonではなくなりますが。
from PIL import ImageGrab
class Window(Tkinter.Tk):
def saveImage(self):
filepath = 'c://Temp/dem/capture%05d.png' % self.dem.step
img = ImageGrab.grab()
s,x,y = self.geometry().split('+')
w,h = s.split('x')
w,h,x,y = map(int,[w,h,x,y])
x += 8
y += 30
img = img.crop((x,y,x+w,y+h))
img.save(filepath)