高速化

DEMソルバーをCythonで作成し、高速化します。生Pythonに比べて大体40倍くらい早くなります。

DEMソルバー（Cython）

Cythonのソースコードです。
コンパイル済みのpydファイル（64bit）をダウンロードしたい場合は、コチラ。

# -*- coding: utf-8 -*-
 
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 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 int_no[9][9]
cdef int area[4] # 解析範囲
cdef int pe_count #粒子の数
cdef int le_count #線要素の数
cdef double dt = 0.001 #計算間隔
cdef int st = 0 #ステップ
 
cdef double rho = 10 #粒子間密度
 
cdef void next_step():
    global st
    reset_force()
    calc_force()
    update_coord()
    st += 1
 
cdef void reset_force():
    cdef int i
    for i in range(pe_count):
        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 = int_no[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 calc_force():
    #2粒子間の接触判定
    cdef double lx,ly,ld,cos_a,sin_a
    cdef int i,j,n
    for i in range(pe_count):
        for j in range(i+1, pe_count):
            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
                calc_force_par2par(i,j,cos_a,sin_a)
            else:
                pe[i].en[j] = 0.0
                pe[i].es[j] = 0.0
 
    #粒子と線の接触判定
    cdef bool hit
    cdef double x,y,a,d,b,s
    for i in range(pe_count):
        for j in range(le_count):
            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)
                s = sqrt(s**2 + pe[i].r**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
                calc_force_line2par(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(pe_count):
        pe[i].fy += -G*pe[i].m #重力
 
cdef void calc_force_par2par(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)
 
    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 calc_force_line2par(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 update_coord():
    cdef double ax,ay,aa
    cdef int i
    for i in range(pe_count):
        #位置更新（オイラー差分）
        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 set_delta_time(sec):
    global dt
    dt = sec
 
def set_number_of_particle(n):
    global pe_count,pe
    pe_count = n
    pe = <Particle*> PyMem_Malloc(n * sizeof(Particle))
    if not pe:
        raise MemoryError()
 
def set_particle(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 set_number_of_line(n):
    global le_count,le
    le_count = n
    le = <Line*> PyMem_Malloc(n * sizeof(Line))
    if not le:
        raise MemoryError()
 
def set_line(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
 
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_obj
 
def set_area(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 set_interface(mat1,mat2,inf_no,inf_obj):
    int_no[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():
    global pe,le
    cdef int i,j,n
    n = pe_count + le_count + 4
    for i in range(pe_count):
        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 = <double*> PyMem_Malloc(n * sizeof(double))
        pe[i].es = <double*> PyMem_Malloc(n * sizeof(double))
        for j in range(n):
            pe[i].en[j] = 0
            pe[i].es[j] = 0
    for i in range(le_count):
        le[i].etype = 2
        le[i].n = pe_count+i
 
def setup():
    global pe,le
    cdef int i
    #粒子要素
    for i in range(pe_count):
        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 calc_step(int n=1):
    cdef int i
    for i in range(n):
        next_step()

ユーザーインターフェイス（Python）

Cythonで作成したDEMソルバーを動かすPythonのソースコードです。
実行すると、100個以上のの球がコロコロ転がる描写が実行させるはずです。

# -*- coding: utf-8 -*-
 
print('読み込み中...',end='')
 
import sys
import math
import random
import time
import tkinter
 
import dem
 
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=0,y1=0,x2=0,y2=0):
        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):
        # setup
        self.area = [5,295,5,195]
 
        # lines
        lines = []
        lines.append(Line(100,100,300,150))
        lines.append(Line(10,80,160,50))
        a = self.area
        lines.append(Line(a[0],a[2],a[1],a[2]))
        lines.append(Line(a[0],a[3],a[1],a[3]))
        lines.append(Line(a[0],a[2],a[0],a[3]))
        lines.append(Line(a[1],a[2],a[1],a[3]))
        self.line_count = len(lines)
 
        # particle
        pars = []
        for x in range(40,290,2):
            for y in range(145,190,2):
                if self._is_touch_to_particles(x,y,5,pars):
                    continue
                if self._is_touch_to_lines(x,y,5,lines):
                    continue
                pars.append(Particle(x,y))
        self.par_count = 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.set_delta_time(0.01)
        dem.set_area(*self.area)
        dem.set_number_of_particle(self.par_count)
        dem.set_number_of_line(len(lines))
        dem.initialize()
        for i,l in enumerate(lines):
            dem.set_line(i,l)
        for i,p in enumerate(pars):
                dem.set_particle(i,p)
        dem.set_interface(1,1,0,inf[1][1])
        dem.set_interface(1,2,1,inf[1][2])
        dem.setup()
 
    def _is_touch_to_particles(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 _is_touch_to_lines(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.par_count):
            p = dem.particle(i,Particle())
            pars.append(p)
        return pars
 
    def lines(self):
        lines = []
        for i in range(self.line_count):
            l = dem.line(i,Line())
            lines.append(l)
        return lines
 
 
class Window(tkinter.Tk):
 
    def __init__(self):
        self.loop_n = 0
 
        print('初期設定中...',end='')
        super().__init__()
        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()
        for l in self.dem_ui.lines():
            xy = self.view_coord([l.x1,l.y1,l.x2,l.y2])
            self.canvas.create_line(xy,width=1)
 
        self.redraw()
        self.update_idletasks()
 
        print('完了')
        print('粒子要素数： %d ' % self.dem_ui.par_count)
        print('解析開始')
 
    def calcloop(self):
        if self.loop_n % 10 == 0:
            print('Step %d' % dem.step())
        if self.loop_n % 1 == 0:
            self.redraw()
        if dem.step() >= 5000:
            print('解析終了.設定最大ループに達しました')
        else:
            self.after(0,self.calcloop)
        dem.calc_step(10)
        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.view_coord([p.x-p.r,p.y-p.r])
            x2,y2 = self.view_coord([p.x+p.r,p.y+p.r])
            self.canvas.create_oval(x1,y1,x2,y2,tags='elem')
 
            x1,y1 = self.view_coord([p.x,p.y])
            x2,y2 = self.view_coord([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 view_coord(self,coords,offset=(0,0)):
        s = 1.0 # 表示倍率
        h = 200 #表示画面高さ
        w = 300 #表示画面幅
        x_offset = 0
        y_offset = 0
        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 main():
    w = Window()
    w.after(0,w.calcloop)
    w.mainloop()
 
print('完了') #読み込み
 
if __name__ == '__main__':
    main()