3Dグラフィック、プログラミング、ものづくり・・・、創作活動を楽しもう!

DEMの勉強3

3といっても新しいことはしていないのだけど。

前回記事の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:             #法線力がなければ、せん断力は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
        p2.fx +=  hn*cos_a - hs*sin_a
        p2.fy +=  hn*sin_a + hs*cos_a
        p2.fm -=  p2.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)

Updated: 2018年1月3日 — 21:35

コメントを残す

メールアドレスが公開されることはありません。

TibLab Blog © 2014 Frontier Theme