DEMの勉強8

いよいよモデルいじり。粒子を土のモデルにしていきます。
モデルは色々あれど、一番最初に検索したら出てきたこの論文を参考にすることにします。

 

地盤材料の破壊基準を表現するためのシンプルな個別要素モデル

クリックして71.pdfにアクセス

 

ポイントは2つ。モールクーロンの破壊基準を適用することと、転がり抵抗を導入すること。この内、転がり抵抗を導入してみました。

 

転がり抵抗は、実際の粒子は角ばっていて、コロコロ転がる事はないのですが、円形要素では、コロコロ転がっていってしまうので、それを防ぐために導入します。要は、円形要素で角ばった砂の粒子を模擬するためのものです。

さて、ちゃんと組めているか?


(左:転がり抵抗無し 右:転がり抵抗有り)

 

一瞬、戻っているような・・・気がするけど、コードがおかしいのか、せん断バネの反発で変に見えるだけか。とりあえず、抵抗力は発揮されているみたい。

さあ、沢山粒子のあるパターンで、安息角に違いが出るか確かめたい。

逆に、流れ出す粒子の数が増えてしまったなり。はて。イメージと違う。

 

ちなみに完全に回転を止めてみると・・・、

いくつも角が生えて、非現実的な面白い感じになります。

 

さて、どうやってちゃんとできているか確かめよう・・。

dem.pyx (cython)

# -*- coding: utf-8 -*-

from libc.math cimport sqrt,sin,cos,tan,atan2,fabs,ceil,floor
from cpython cimport bool
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free

cdef double PI = 3.1415926535 #円周率
cdef double G = 9.80665 #重力加速度

class _Particle:
    pass
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 #慣性モーメント
    double ph #せん断抵抗角
    double t #粘着力パラメータ(N/m)
    int *nearList #近傍リスト
    int nearCount #近傍リスト数
    int *co #ボンド接着

class _Line:
    pass
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

class _Interface:
    pass
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 elCount #全要素数
cdef int peCount #粒子の数
cdef int leCount #線要素の数
cdef double dt = 0.001 #計算間隔
cdef int st = 0 #ステップ

cdef double rho = 10 #粒子間密度

def _nextStep():
    pass
cdef void nextStep():
    global st,nearUpdate
    if nearUpdate == True:
        cellRegist()
    resetForce()
    calcForce()
    if nearUpdate == True:
        nearUpdate = False
    updateCoord()
    st += 1
        
def _resetForce():
    pass
cdef void resetForce():
    cdef int i
    for i in range(peCount):
        pe[i].fx = 0
        pe[i].fy = 0
        pe[i].fm = 0

def _interface():
    pass
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

def _calcForce():
    pass
cdef void calcForce():
    #2粒子間の接触判定
    cdef double lx,ly,ld,cos_a,sin_a
    cdef int i,j,n,m
    for i in range(peCount):
        #粒子登録法による近傍リストの絞込み
        if nearUpdate == True:
            m = 0
            for n in range(pe[i].nearCount):
                j = pe[i].nearList[n]
                rc = pe[i].r + pe[j].r + nearAlpha * r_min*2
                lx = pe[j].x - pe[i].x
                ly = pe[j].y - pe[i].y
                ld = (lx**2+ly**2)**0.5
                if ld < rc:
                    nearList[m] = j
                    m += 1
            pe[i].nearCount = m
            for n in range(pe[i].nearCount):
                pe[i].nearList[n] = nearList[n]
        
        #for j in range(peCount):
        for n in range(pe[i].nearCount):
            j = pe[i].nearList[n]
            if i == j:
                continue
            lx = pe[j].x - pe[i].x
            ly = pe[j].y - pe[i].y
            ld = (lx**2+ly**2)**0.5
                 
            if ld<(pe[i].r+pe[j].r):
                cos_a = lx/ld
                sin_a = ly/ld
                force2par(i,j,cos_a,sin_a,ld*pe[i].r/(pe[i].r+pe[j].r))
            else:
                pe[i].en[j] = 0.0
                pe[i].es[j] = 0.0
        
    #粒子と壁の接触判定
    cdef double cond[4]
    cdef double xs[4]
    cdef double 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
                force2par(i,n,cos_a,sin_a,ld)
            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
                force2par(i,le[j].n,cos_a,sin_a,ld)
            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 #重力

def _force2par():
    pass
cdef void force2par(int i,int j,double cos_a,double sin_a,double rd):
    cdef int m,n
    cdef int mat[2]
    cdef double dx[2]
    cdef double dy[2]
    cdef double da[2]
    cdef double r[2]
    cdef double un,us,vn,vs,hn,hs
    cdef Interface inf
    
    for m in range(2):
        if m == 0 : n = i
        else:       n = j
        if n < peCount:
            dx[m] = pe[n].dx
            dy[m] = pe[n].dy
            da[m] = pe[n].da
            r[m] = pe[n].r
            mat[m] = pe[n].etype
        else:
            dx[m] = 0
            dy[m] = 0
            da[m] = 0
            r[m] = 0
            mat[m] = 2
                
    #相対的変位増分
    un = +(dx[0]-dx[1])*cos_a+(dy[0]-dy[1])*sin_a
    us = -(dx[0]-dx[1])*sin_a+(dy[0]-dy[1])*cos_a+(r[0]*da[0]+r[1]*da[1])
    #相対的速度増分
    #vn = +(pe[0].vx-pe[1].vx)*cos_a+(pe[0].vy-pe[1].vy)*sin_a
    #vs = -(pe[0].vx-pe[1].vx)*sin_a+(pe[0].vy-pe[1].vy)*cos_a+(pe[0].r*pe[0].va+pe[1].r*pe[1].va)
    
    inf = interface(mat[0],mat[1])
    
    #合力(局所座標系)
    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 #時計回りが正
    
    # 粒子間ボンドモデル
    cdef double bs,tu
    if False:# i < peCount and j < peCount:
        tu = pe[i].t * (pe[i].r + pe[j].r)
        if -tu < hn: pe[i].co[j] = False print('connect break',i,j) if pe[i].co[j] == True: bs = (hn + tu) * tan(pe[i].ph) else: bs = hn * tan(pe[i].ph) if hs > 0:
            if hs > bs:
                hs = bs
        else: # hs < 0
            if hs < bs:
                hs = bs
    else:
        if hn <= 0.0: #法線力がなければ、せん断力は0 hs = 0.0 elif fabs(hs) >= inf.frc*hn:
            #摩擦力以上のせん断力は働かない
            hs = inf.frc*fabs(hn)*hs/fabs(hs)     
    
    #転がり摩擦抵抗
    cdef double mr
    cdef double b = 0.3
    mr = hn * b * rd
    if fabs(hs)*pe[i].r < fabs(mr): #逆回転するようなMrは発生しない mr = hs*pe[i].r else: #打ち消す方向に作用 if hs > 0: mr =  fabs(mr)
        else:      mr = -fabs(mr)
    
    #全体合力(全体座標系)
    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 - mr

def _updateCoord():
    pass
cdef void updateCoord():
    global xbook
    cdef double ax,ay,aa,v,vmax
    vmax = 0
    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
        #以下近傍リスト用処理
        v = sqrt(pe[i].vx**2+pe[i].vy**2)
        if vmax < v: vmax = v xbook += vmax * dt if xbook > nearAlpha*r_min:
        global nearUpdate
        nearUpdate = True
        xbook = 0

#セル分割法と粒子登録法
##cell = []
##cellCount = []
cdef int ***cell
cdef int **cellCount
cdef int cellInfo[3]
cdef bool nearUpdate = False
cdef double nearAlpha = 0.1
cdef double xbook
cdef int *nearList 

def _cellInit():
    pass
cdef void cellInit():
    global cell,cellCount,nearList
    cell_width = r_max*2 + nearAlpha * r_min*2
    xn = int(ceil((area[1]-area[0])/cell_width))
    yn = int(ceil((area[3]-area[2])/cell_width))
    cn = int(ceil(cell_width/r_min)+1)**2 #+1予備(重なりを考慮)
##    cell = [[[-1 for i in range(cn)] for j in range(yn)] for k in range(xn)]
##    cellCount = [[0 for i in range(yn)] for j in range(xn)]
    cell = <int***> PyMem_Malloc(xn * sizeof(int*))
    cellCount = <int**> PyMem_Malloc(xn * sizeof(int*))
    for i in range(xn):
        cell[i] = <int**> PyMem_Malloc(yn * sizeof(int*))
        cellCount[i] = <int*> PyMem_Malloc(yn * sizeof(int))
        for j in range(yn):
            cell[i][j] = <int*> PyMem_Malloc(cn * sizeof(int))
    cellInfo[0] = xn #列数
    cellInfo[1] = yn #行数
    cellInfo[2] = cn #セル格納最大数
    cellReset()
    nearList = <int*> PyMem_Malloc(cn * sizeof(int))
    
def _cellRegist():
    pass
cdef void cellRegist():
    cdef double cell_width
    cdef int xn,yn,cn,n,i,k,j,l
    cellReset()
    cell_width = r_max*2 + nearAlpha * r_min*2
    #セルに要素No格納
    for i in range(peCount):
        xn = int(floor(pe[i].x/cell_width))
        yn = int(floor(pe[i].y/cell_width))
        if xn<0 or cellInfo[0]<=xn or yn<0 or cellInfo[1]<=yn:
            #print('Error! particle coord is out of range.')
            #print('Particle n=%d x=%0.3f y=%0.3f' % (i,pe[i].x,pe[i].y))
            continue 
        cn = cellCount[xn][yn]
        cell[xn][yn][cn] = i
        cellCount[xn][yn] += 1
    #近傍リスト作成
    for i in range(peCount):
        xn = int(floor(pe[i].x/cell_width))
        yn = int(floor(pe[i].y/cell_width))
        n = 0
        for j in range(xn-1,xn+2):
            for k in range(yn-1,yn+2):
                if j < 0 or cellInfo[0] <= j:
                    continue
                if k < 0 or cellInfo[1] <= k:
                    continue
                for l in range(cellCount[j][k]):
                    pe[i].nearList[n] = cell[j][k][l]
                    n += 1
        pe[i].nearCount = n
        
cdef void cellReset():
    cdef int i,j,k
    for i in range(cellInfo[0]):
        for j in range(cellInfo[1]):
            cellCount[i][j] = 0
            for k in range(cellInfo[2]):
                cell[i][j][k] = -1

# -------------------------
# Pythonからの設定用
# -------------------------

def setNumberOfParticle(n):
    global peCount,pe
    peCount = n
    pe = <Particle*> PyMem_Malloc(n * sizeof(Particle))
    if not pe:
        raise MemoryError()

def setNumberOfLine(n):
    global leCount,le
    leCount = n
    le = <Line*> PyMem_Malloc(n * sizeof(Line))
    if not le:
        raise MemoryError()

def initialize():
    global elCount
    elCount = peCount + leCount + 4
    cdef int i,j
    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 = <double*> PyMem_Malloc(elCount * sizeof(double))
        pe[i].es = <double*> PyMem_Malloc(elCount * sizeof(double))
        pe[i].co = <int*> PyMem_Malloc(elCount * sizeof(int))
        for j in range(elCount):
            pe[i].en[j] = 0
            pe[i].es[j] = 0
            pe[i].co[j] = True
    for i in range(leCount):
        le[i].etype = 2
        le[i].n = peCount+i

def setDeltaTime(sec):
    global dt
    dt = sec

def setParticle(pe_no,pe_obj):
    pe[pe_no].r = pe_obj.r
    pe[pe_no].x = pe_obj.x
    pe[pe_no].y = pe_obj.y
    pe[pe_no].a = pe_obj.a
    pe[pe_no].dx = pe_obj.dx
    pe[pe_no].dy = pe_obj.dy
    pe[pe_no].da = pe_obj.da
    pe[pe_no].vx = pe_obj.vx
    pe[pe_no].vy = pe_obj.vy
    pe[pe_no].va = pe_obj.va 
    pe[pe_no].fx = pe_obj.fx 
    pe[pe_no].fy = pe_obj.fy
    pe[pe_no].fm = pe_obj.fm
    pe[pe_no].m = pe_obj.m
    pe[pe_no].Ir = pe_obj.Ir
    #pe[pe_no].ph = pe_obj.ph
    #pe[pe_no].t = pe_obj.t    
    cdef int i
    for i in range(elCount):
        if i < len(pe_obj.en): pe[pe_no].en[i] = pe_obj.en[i] pe[pe_no].es[i] = pe_obj.es[i] def particle(pe_no,pe_obj): pe_obj.n = pe[pe_no].n pe_obj.r = pe[pe_no].r pe_obj.x = pe[pe_no].x pe_obj.y = pe[pe_no].y pe_obj.a = pe[pe_no].a pe_obj.dx = pe[pe_no].dx pe_obj.dy = pe[pe_no].dy pe_obj.da = pe[pe_no].da pe_obj.vx = pe[pe_no].vx pe_obj.vy = pe[pe_no].vy pe_obj.va = pe[pe_no].va pe_obj.fx = pe[pe_no].fx pe_obj.fy = pe[pe_no].fy pe_obj.fm = pe[pe_no].fm pe_obj.m = pe[pe_no].m pe_obj.Ir = pe[pe_no].Ir pe_obj.ph = pe[pe_no].ph pe_obj.t = pe[pe_no].t pe_obj.en = [0.0 for i in range(elCount)] pe_obj.es = [0.0 for i in range(elCount)] for i in range(elCount): pe_obj.en[i] = pe[pe_no].en[i] pe_obj.es[i] = pe[pe_no].es[i] return pe_obj 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 cdef double r_max #粒子の最大半径 cdef double r_min #粒子の最小半径 def setup(): global r_max,r_min cdef int i r_max = 0 r_min = float('inf') for i in range(peCount): if pe[i].r > r_max:
            r_max = pe[i].r
        if pe[i].r < r_min:
            r_min = pe[i].r
    cellInit()
    pn = cellInfo[2] * 9
    for i in range(peCount):
        pe[i].nearList = <int*> PyMem_Malloc(pn * sizeof(int))
    cellRegist()

def step():
    return st

def calcStep(int n=1):
    cdef int i
    for i in range(n):
        nextStep()

 

dem_ui.py

# -*- coding: utf-8 -*-

print u'読み込み中...',

import sys
import math
import random
import time
import cPickle
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(object):
    def __init__(self,x=0,y=0,vx=0,vy=0):
        #要素共通
        self.etype = 1 #要素タイプ
        self.n = 0 #要素No.
        self.r = 5.0 #半径
        self.x = x #X座標
        self.y = y #Y座標
        self.a = 0.0 #角度
        self.dx = 0.0 #X方向増加量
        self.dy = 0.0 #Y方向増加量
        self.da = 0.0 #角度増加量
        self.vx = vx #X方向速度
        self.vy = vy #Y方向速度
        self.va = 0.0 #角速度
        self.fx = 0.0 #X方向節点力
        self.fy = 0.0 #Y方向節点力
        self.fm = 0.0 #回転モーメント
        self.en = [] #弾性力(直方向)
        self.es = [] #弾性力(せん断方向)
    
        #粒子専用
        rho = 10
        self.m = 4.0/3.0*math.pi*rho*self.r**3 # 質量
        self.Ir = math.pi*rho*self.r**4/2.0 #慣性モーメント
        #土粒子
        self.ph = math.radians(30) #せん断抵抗角
        self.t = 20 #粘着力を決めるパラメータ N/m

class Line(object):
    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._lines = []
        
        self._config()
        #self._test() 
        
    def _test(self):
        self.area = [5,295,5,195]
        self.parCount = 2
        dem.setArea(*self.area)
        dem.setDeltaTime(0.01)
        dem.setNumberOfParticle(self.parCount)
        dem.setNumberOfLine(1)
        dem.initialize()
        p1 = Particle(100,60)
        p1.r = 10
        p2 = Particle(100,90)
        p2.r = 10
        dem.setParticle(0,p1)
        dem.setParticle(1,p2)
        l1 = Line(0,30,200,10)
        self._lines.append(l1)
        dem.setLine(0,l1)
        #粒子同士
        inf = [[Interface() for i in range(9)] for j in range(9)]
        inf[1][1].kn  = 100000 #弾性係数(法線方向)
        inf[1][1].etan= 50000 #粘性係数(法線方向)
        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 = 100
        dem.setInterface(1,1,0,inf[1][1])
        dem.setInterface(1,2,1,inf[1][2])
        dem.setup()
        
    def _config(self):
        
        dt = 0.01
        #area = [5,295,5,195]
        area = [5,295,-500,195]
        step = 2
        
        # lines
        lines = []
        lines.append(Line(5,40,150,50))
        if step <= 1: lines.append(Line(150,50,150,195)) else: lines.append(Line(150,50,150,60)) # particle pars = [] if step == 0: _area = [10,130,45,190] pars = self._setParInArea(_area,pars,lines) elif step == 1: pars = self.load() #さらに粒子追加 #_area = [150,290,10,90] _area = [10,140,130,190] add_pars = self._setParInArea(_area,pars,lines) #弾性力調整 ad = [0 for i in range(len(add_pars))] for p in pars: p.en = p.en[:len(pars)] + ad + p.en[len(pars):] p.es = p.es[:len(pars)] + ad + p.es[len(pars):] pars = pars + add_pars else: pars = self.load() # interface inf1 = Interface() #粒子同士 inf1.kn = 100000 #弾性係数(法線方向) inf1.etan= 50000 #粘性係数(法線方向) inf1.ks = 5000 #弾性係数(せん断方向) inf1.etas= 1000 #粘性係数(せん断方向) inf1.frc = 10 #摩擦係数 #粒子と線要素 inf2 = Interface() inf2.kn = 500000 inf2.etan= 10000 inf2.ks = 1000 inf2.etas= 900 inf2.frc = 100 infs = [] infs.append([1,1,0,inf1]) infs.append([1,2,1,inf2]) self.setup(dt,area,pars,lines,infs) def _setParInArea(self,area,pars,lines,pr=5): # area = [x_min,x_max,y_min,y_max] add_pars = [] for x in range(area[0],area[1],1): for y in range(area[2],area[3],1): _pars = pars + add_pars if self._hitParticle(x,y,pr,_pars): continue if self._hitLine(x,y,pr,lines): continue add_pars.append(Particle(x,y)) return add_pars 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 setup(self,dt,area,pars,lines,infs): #setup dem dem.setDeltaTime(dt) dem.setArea(*area) dem.setNumberOfParticle(len(pars)) dem.setNumberOfLine(len(lines)) dem.initialize() for i,l in enumerate(lines): dem.setLine(i,l) for i,p in enumerate(pars): dem.setParticle(i,p) for v in infs: dem.setInterface(*v) dem.setup() self.area = area self.parCount = len(pars) self._lines = lines 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 def save(self,fp='init.dem'): f = open(fp,'wb') f.write(cPickle.dumps(self.particles())) f.close() print('save particles data at '+fp) def load(self,fp='init.dem'): f = open(fp,'rb') pars = cPickle.loads(f.read()) f.close() return pars class Window(Tkinter.Tk): def __init__(self): print u'初期設定中...', self.loop_n = 0 self.width = 300 self.height = 200 self.scale = 1.0 Tkinter.Tk.__init__(self) self.canvas = Tkinter.Canvas(self, bg="white") self.canvas.pack(fill=Tkinter.BOTH,expand=True) self.geometry('%dx%d' % (self.width,self.height)) self.title('DEM') self.dem_ui = DEM_UI() a = self.dem_ui.area area = [a[0],a[2],a[1],a[2],a[1],a[3],a[0],a[3],a[0],a[2]] xy = self.viewCoord(area) self.canvas.create_line(xy) 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'完了') print(u'粒子要素数: %d ' % self.dem_ui.parCount) print(u'解析開始') def calcloop(self): if self.loop_n == 1: self.saveCalcTime('start') pass 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 >= 500:
            print(u'解析終了.設定最大ループに達しました')
            self.saveCalcTime('finish')
            #self.dem_ui.save()
        else:
            self.after(0,self.calcloop)
        dem.calcStep(10)
        self.loop_n += 1
        self.update_idletasks()
 
    def redraw(self):
        self.canvas.delete('elem')
        h = self.height
        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 = self.scale # 表示倍率
        h = self.height #表示画面高さ
        w = self.width #表示画面幅
        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()
    
def test():
    dem_ui = DEM_UI()
    en1 = dem_ui.particles()[100].es
    pars = dem_ui.load()
    en2 = pars[100].es
    for i in range(len(en1)):
        v = en1[i] - en2[i]
        #if abs(v) > 1:
        print(i,v)#,en1[i],en2[i])

print u'完了'

if __name__ == '__main__':
    main()

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です