前回のプログラムを、バク直したり、使いたくないかったポインタ変数使ってみたりして、一区切りついた感じなので、投稿しときます。
この先は、土木屋っぽく土を疑似するように必要な機能を実装していく予定。
Cython:dem.pyx [dem.pyd]
# -*- coding: utf-8 -*-
#Cython仕様で一から作ってみる
from libc.math cimport sqrt,sin,cos,atan2,fabs
from cpython cimport bool
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
cdef double PI = 3.1415926535 #円周率
cdef double G = 9.80665 #重力加速度
cdef struct Particle:
#要素共通
int etype #要素タイプ
int n #要素No.
double r#半径
double x#X座標
double y#Y座標
double a#角度
double dx#X方向増加量
double dy#Y方向増加量
double da#角度増加量
double vx#X方向速度
double vy#Y方向速度
double va#角速度
double fy
double fx
double fm
double *en #弾性力（直方向）
double *es #弾性力（せん断方向）
#粒子専用
double m# 質量
double Ir#慣性モーメント
cdef struct Line:
#要素共通
int etype #要素タイプ
int n #要素No.
double r#半径
double x#X座標
double y#Y座標
double a#角度
double dx#X方向増加量
double dy#Y方向増加量
double da#角度増加量
double vx#X方向速度
double vy#Y方向速度
double va#角速度
double fy
double fx
double fm
#double *en #弾性力（直方向）
#double *es #弾性力（せん断方向）
#線要素専用
double x1
double y1
double x2
double y2
cdef struct Interface:
double kn #弾性係数（法線方向）
double etan#粘性係数（法線方向）
double ks#弾性係数（せん断方向）
double etas#弾性係数（せん断方向）
double frc#摩擦係数
#グローバル変数
cdef Particle *pe
cdef Line *le
cdef Interface infs[9]
cdef int infNo[9][9]
cdef int area[4] # 解析範囲
cdef int peCount #粒子の数
cdef int leCount #線要素の数
cdef double dt = 0.001 #計算間隔
cdef int st = 0 #ステップ
cdef double rho = 10 #粒子間密度
cdef void nextStep():
resetForce()
calcForce()
updateCoord()
global st
st += 1
cdef void resetForce():
cdef int i
for i in range(peCount):
pe[i].fx = 0
pe[i].fy = 0
pe[i].fm = 0
cdef Interface interface(int etype1,int etype2):
cdef Interface inf
cdef int n
n = infNo[etype1][etype2]
inf.kn = infs[n].kn
inf.etan = infs[n].etan
inf.ks = infs[n].ks
inf.etas = infs[n].etas
inf.frc = infs[n].frc
return inf
cdef void calcForce():
#2粒子間の接触判定
cdef double lx,ly,ld,cos_a,sin_a
cdef int i,j,n
for i in range(peCount):
for j in range(i+1,peCount):
lx = pe[j].x - pe[i].x
ly = pe[j].y - pe[i].y
ld = (lx**2+ly**2)**0.5
if (pe[i].r+pe[j].r)>ld:
cos_a = lx/ld
sin_a = ly/ld
forcePar2par(i,j,cos_a,sin_a)
else:
pe[i].en[j] = 0.0
pe[i].es[j] = 0.0
#粒子と壁の接触判定
cdef double cond[4],xs[4],ys[4]
for i in range(peCount):
cond = [pe[i].x-pe[i].r-area[0],-(pe[i].x+pe[i].r-area[1]),
pe[i].y-pe[i].r-area[2],-(pe[i].y+pe[i].r-area[3])]
xs = [area[0],area[1],pe[i].x,pe[i].x]
ys = [pe[i].y,pe[i].y,area[2],area[3]]
for j in range(4):
n = peCount + leCount + j
if cond[j] < 0:
lx = xs[j] - pe[i].x
ly = ys[j] - pe[i].y
ld = sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
forceLine2par(i,n,2,cos_a,sin_a)
else:
pe[i].en[n] = 0.0
pe[i].es[n] = 0.0
#粒子と線の接触判定
cdef bool hit
cdef double x,y,a,d,b,s
for i in range(peCount):
for j in range(leCount):
hit = False
th0 = atan2(le[j].y2-le[j].y1, le[j].x2-le[j].x1)
th1 = atan2(pe[i].y -le[j].y1, pe[i].x -le[j].x1)
a = sqrt((pe[i].x-le[j].x1)**2+(pe[i].y-le[j].y1)**2)
d = fabs(a*sin(th1-th0))
if d < pe[i].r:
b = sqrt((pe[i].x -le[j].x2)**2+(pe[i].y -le[j].y2)**2)
s = sqrt((le[j].x2-le[j].x1)**2+(le[j].y2-le[j].y1)**2)
if a < s and b < s:
s = sqrt(a**2-d**2)
x = le[j].x1 + s*cos(th0)
y = le[j].y1 + s*sin(th0)
hit = True
elif a < b and a < pe[i].r:
x = le[j].x1
y = le[j].y1
hit = True
elif b < pe[i].r:
x = le[j].x2
y = le[j].y2
hit = True
if hit:
lx = x - pe[i].x
ly = y - pe[i].y
ld = sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
forceLine2par(i,le[j].n,2,cos_a,sin_a)
else:
pe[i].en[le[j].n] = 0.0
pe[i].es[le[j].n] = 0.0
#外力
for i in range(peCount):
pe[i].fy += -G*pe[i].m #重力
cdef void forcePar2par(int i,int j,double cos_a,double sin_a):
cdef double un,us,vn,vs,hn,hs
cdef Interface inf
#相対的変位増分
un = +(pe[i].dx-pe[j].dx)*cos_a+(pe[i].dy-pe[j].dy)*sin_a
us = -(pe[i].dx-pe[j].dx)*sin_a+(pe[i].dy-pe[j].dy)*cos_a+(pe[i].r*pe[i].da+pe[j].r*pe[j].da)
#相対的速度増分
#vn = +(pe[i].vx-pe[j].vx)*cos_a+(pe[i].vy-pe[j].vy)*sin_a
#vs = -(pe[i].vx-pe[j].vx)*sin_a+(pe[i].vy-pe[j].vy)*cos_a+(pe[i].r*pe[i].va+pe[j].r*pe[j].va)
inf = interface(pe[i].etype,pe[j].etype)
#合力（局所座標系）
pe[i].en[j] += inf.kn*un
pe[i].es[j] += inf.ks*us
hn = pe[i].en[j] + inf.etan*un/dt
hs = pe[i].es[j] + inf.etas*us/dt
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:
#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn)*hs/fabs(hs)
#全体合力（全体座標系）
pe[i].fx += -hn*cos_a + hs*sin_a
pe[i].fy += -hn*sin_a - hs*cos_a
pe[i].fm -= pe[i].r*hs
pe[j].fx += hn*cos_a - hs*sin_a
pe[j].fy += hn*sin_a + hs*cos_a
pe[j].fm -= pe[j].r*hs
cdef void forceLine2par(int i,int ln,int mat,double cos_a, double sin_a):
cdef double un,us,vn,vs,hn,hs
cdef Interface inf
#相対的変位増分
un = +pe[i].dx*cos_a+pe[i].dy*sin_a
us = -pe[i].dx*sin_a+pe[i].dy*cos_a+pe[i].r*pe[i].da
#相対的速度増分
vn = +pe[i].vx*cos_a+pe[i].vy*sin_a
vs = -pe[i].vx*sin_a+pe[i].vy*cos_a+pe[i].r*pe[i].va
inf = interface(pe[i].etype,mat)
#合力（局所座標系）
pe[i].en[ln] += inf.kn*un
pe[i].es[ln] += inf.ks*us
hn = pe[i].en[ln] + inf.etan*vn
hs = pe[i].es[ln] + inf.etas*vs
if hn <= 0.0:#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn)*hs/fabs(hs)
#全体合力（全体座標系）
pe[i].fx += -hn*cos_a + hs*sin_a
pe[i].fy += -hn*sin_a - hs*cos_a
pe[i].fm -= pe[i].r*hs
cdef void updateCoord():
cdef double ax,ay,aa
cdef int i
for i in range(peCount):
#位置更新（オイラー差分）
ax = pe[i].fx/pe[i].m
ay = pe[i].fy/pe[i].m
aa = pe[i].fm/pe[i].Ir
pe[i].vx += ax*dt
pe[i].vy += ay*dt
pe[i].va += aa*dt
pe[i].dx = pe[i].vx*dt
pe[i].dy = pe[i].vy*dt
pe[i].da = pe[i].va*dt
pe[i].x += pe[i].dx
pe[i].y += pe[i].dy
pe[i].a += pe[i].da
# -------------------------
# Pythonからの設定用
# -------------------------
def setDeltaTime(sec):
global dt
dt = sec
def setNumberOfParticle(n):
global peCount,pe
peCount = n
pe = PyMem_Malloc(n * sizeof(Particle))
if not pe:
raise MemoryError()
def setParticle(pe_no,pe_obj):
pe[pe_no].x = pe_obj.x
pe[pe_no].y = pe_obj.y
def particle(pe_no,pe_obj):
pe_obj.x = pe[pe_no].x
pe_obj.y = pe[pe_no].y
pe_obj.a = pe[pe_no].a
return pe_obj
def setNumberOfLine(n):
global leCount,le
leCount = n
le = PyMem_Malloc(n * sizeof(Line))
if not le:
raise MemoryError()
def setLine(l_no,l_obj):
le[l_no].x1 = l_obj.x1
le[l_no].y1 = l_obj.y1
le[l_no].x2 = l_obj.x2
le[l_no].y2 = l_obj.y2
return l_obj
def line(l_no,l_obj):
l_obj.x1 = le[l_no].x1
l_obj.y1 = le[l_no].y1
l_obj.x2 = le[l_no].x2
l_obj.y2 = le[l_no].y2
return l_no
def setArea(x_min,x_max,y_min,y_max):
global area
area[0] = x_min
area[1] = x_max
area[2] = y_min
area[3] = y_max
def setInterface(mat1,mat2,inf_no,inf_obj):
infNo[mat1][mat2] = inf_no
infs[inf_no].kn = inf_obj.kn
infs[inf_no].etan = inf_obj.etan
infs[inf_no].ks = inf_obj.ks
infs[inf_no].etas = inf_obj.etas
infs[inf_no].frc = inf_obj.frc
def initialize():
cdef int i,j,n
n = peCount + leCount + 4
for i in range(peCount):
pe[i].etype = 1
pe[i].n = i
pe[i].x = 0
pe[i].y = 0
pe[i].r = 0
pe[i].a = 0
pe[i].dx = 0
pe[i].dy = 0
pe[i].da = 0
pe[i].vx = 0
pe[i].vy = 0
pe[i].va = 0
pe[i].fx = 0
pe[i].fy = 0
pe[i].fm = 0
pe[i].m = 0
pe[i].Ir = 0
pe[i].en = PyMem_Malloc(n * sizeof(double))
pe[i].es = PyMem_Malloc(n * sizeof(double))
for j in range(n):
pe[i].en[j] = 0
pe[i].es[j] = 0
for i in range(leCount):
le[i].etype = 2
le[i].n = peCount+i
def setup():
cdef int i
#粒子要素
for i in range(peCount):
pe[i].r = 5.0
pe[i].m = 4.0/3.0*PI*rho*pe[i].r**3 # 質量
pe[i].Ir = PI*rho*pe[i].r**4/2.0 #慣性モーメント
#線要素
def step():
return st
def calcStep(int n=1):
cdef int i
for i in range(n):
nextStep()
Python:dem_ui.py
# -*- coding: utf-8 -*-
print u'読み込み中...',
import sys
import math
import random
import time
import Tkinter
import dem
from PIL import ImageGrab
class Element(object):
def __init__(self):
self.n = 0 #要素No.
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=0,y=0,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方向速度
rho = 10
self.r = 5 #半径
self.m = 4.0/3.0*math.pi*rho*self.r**3 # 質量
self.loop_nr = math.pi*rho*self.r**4/2.0 #慣性モーメント
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(object):
def __init__(self):
self.kn = 0 #弾性係数（法線方向）
self.etan = 0 #粘性係数（法線方向）
self.ks = 0 #弾性係数（せん断方向）
self.etas = 0 #粘性係数（せん断方向）
self.frc = 0 #摩擦係数
class DEM_UI:
def __init__(self):
self._pars = []
self._lines = []
self._setup()
def _setup(self):
self.area = [5,295,5,195]
# lines
self._lines = []
self._lines.append(Line(100,100,300,150))
self._lines.append(Line(10,80,160,50))
# particle
pars = []
for x in range(40,290,2):
for y in range(145,190,2):
if self._hitParticle(x,y,5,pars):
continue
if self._hitLine(x,y,5,self._lines):
continue
pars.append(Particle(x,y))
#if len(pars) >= 1: break
self.parCount = len(pars)
# interface
inf = [[Interface() for i in range(9)] for j in range(9)]
#粒子同士
inf[1][1].kn = 100000 #弾性係数（法線方向）
inf[1][1].etan= 5000 #粘性係数（法線方向）
inf[1][1].ks = 5000 #弾性係数（せん断方向）
inf[1][1].etas= 1000 #粘性係数（せん断方向）
inf[1][1].frc = 10 #摩擦係数
#粒子と線要素
inf[1][2].kn = 500000
inf[1][2].etan= 10000
inf[1][2].ks = 1000
inf[1][2].etas= 900
inf[1][2].frc = 1
#setup dem
dem.setDeltaTime(0.01)
dem.setArea(*self.area)
dem.setNumberOfParticle(self.parCount)
dem.setNumberOfLine(len(self._lines))
dem.initialize()
for i,l in enumerate(self._lines):
dem.setLine(i,l)
for i,p in enumerate(pars):
dem.setParticle(i,p)
dem.setInterface(1,1,0,inf[1][1])
dem.setInterface(1,2,1,inf[1][2])
dem.setup()
print(u'完了')
print(u'粒子要素数： %d ' % self.parCount)
def _hitParticle(self,x,y,r,pars):
hit = False
for p in pars:
lx = p.x - x
ly = p.y - y
ld = (lx**2+ly**2)**0.5
if (p.r+r)>=ld:
hit = True
break
return hit
def _hitLine(self,px,py,pr,lines):
hit = False
for l in lines:
th0 = math.atan2(l.y2-l.y1,l.x2-l.x1)
th1 = math.atan2(py-l.y1,px-l.x1)
a = math.sqrt((px-l.x1)**2+(py-l.y1)**2)
d = abs(a*math.sin(th1-th0))
if d < pr:
b = math.sqrt((px-l.x2)**2+(py-l.y2)**2)
s = math.sqrt((l.x2-l.x1)**2+(l.y2-l.y1)**2)
if a < s and b < s:
hit = True
elif a < b and a < pr:
hit = True
elif b < pr:
hit = True
if hit:
break
return hit
def particles(self):
pars = []
for i in range(self.parCount):
p = dem.particle(i,Particle())
pars.append(p)
return pars
def lines(self):
return self._lines
class Window(Tkinter.Tk):
def __init__(self):
self.loop_n = 0
print u'初期設定中...',
Tkinter.Tk.__init__(self)
self.canvas = Tkinter.Canvas(self, bg="white")
self.canvas.pack(fill=Tkinter.BOTH,expand=True)
self.geometry('300x200')
self.title('DEM')
self.dem_ui = DEM_UI()
a = self.dem_ui.area
self.canvas.create_line(a[0],a[2],a[1],a[2],a[1],a[3],
a[0],a[3],a[0],a[2])
for l in self.dem_ui.lines():
xy = self.viewCoord([l.x1,l.y1,l.x2,l.y2])
self.canvas.create_line(xy,width=1)
self.redraw()
self.update_idletasks()
print(u'解析開始')
def calcloop(self):
dem.calcStep(10)
if self.loop_n == 1:
self.saveCalcTime('start')
if self.loop_n % 5 == 0:
print('Step %d' % dem.step())
if self.loop_n % 1 == 0:
self.redraw()
self.saveImage()
if self.loop_n >= 1000:
self.saveCalcTime('finish')
print(u'解析終了.設定最大ループに達しました')
else:
self.after(0,self.calcloop)
self.update_idletasks()
self.loop_n += 1
def redraw(self):
self.canvas.delete('elem')
h = 200
for p in self.dem_ui.particles():
x1,y1 = self.viewCoord([p.x-p.r,p.y-p.r])
x2,y2 = self.viewCoord([p.x+p.r,p.y+p.r])
self.canvas.create_oval(x1,y1,x2,y2,tags='elem')
x1,y1 = self.viewCoord([p.x,p.y])
x2,y2 = self.viewCoord([p.x+p.r*math.cos(p.a),
p.y+p.r*math.sin(p.a)])
self.canvas.create_line(x1,y1,x2,y2,tags='elem')
def viewCoord(self,coords,offset=(0,0)):
s = 1.0 # 表示倍率
h = 200 #表示画面高さ
w = 300 #表示画面幅
x_offset = 0#int(w/2)
y_offset = 0#int(h/2)
xy_list = []
for i in range(0,len(coords),2):
x = round(s*coords[i])+x_offset
y = round(h-s*coords[i+1])-y_offset
x = x + offset[0]
y = y + offset[1]
xy_list.append(x)
xy_list.append(y)
return xy_list
def saveCalcTime(self,option):
if option == 'start':
self.st_time = time.time()
self.st_step = dem.step()
elif option == 'finish':
now = time.time()
dt = now-self.st_time
ds = dem.step() - self.st_step +1
f = open('calc_time.txt','w')
f.write('START STEP %d\n' % self.st_step)
f.write('START TIME {0}\n'.format(self.st_time))
f.write('END STEP %d\n' % dem.step())
f.write('END TIME {0}\n'.format(now))
f.write('DIFF STEP %d \n' % ds)
f.write('DIFF TIME {0}\n'.format(dt))
f.write('ONE STEP TIME {0}'.format(dt/ds))
f.close()
def saveImage(self):
filepath = 'c://Temp/dem/capture%05d.png' % 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)
def main():
w = Window()
w.after(0,w.calcloop)
w.mainloop()
print u'完了'
if __name__ == '__main__':
main()