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

DEMの勉強2

意味も分からず線とか追加。厳密な当たり判定すらできず、ズルしてみたけど、うまく動いとる。ゲームとかならまだしも、解析でやっちゃいかんでしょー。

# -*- 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:
            #法線力がなければ、せん断力は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:
            #法線力がなければ、せん断力は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
        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(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()

Updated: 2017年10月29日 — 04:49

コメントを残す

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

TibLab Blog © 2014 Frontier Theme