Python コンター図作成

コンター図を作成する必要があったため、そのプログラムを作成した。

作例

プログラムソース

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns


def drawfig(xx,yy,zh,zc):
    q=68.64
    vp1=7.0; dp1=np.sqrt(4*q/np.pi/vp1)
    vp2=2.5; dp2=np.sqrt(4*q/np.pi/vp2)
    vh1=5.0; dh1=np.sqrt(4*q/np.pi/vh1)
    vh2=2.0; dh2=np.sqrt(4*q/np.pi/vh2)
    fsz=16
    xmin=3; xmax=6
    ymin=4; ymax=7
    plt.figure(figsize=(20,10),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'

    plt.subplot(121)
    tstr='Total Head loss Contour of Headrace and Penstock (unit: m)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Penstock diameter (m)')
    plt.ylabel('Headrace diameter (m)')
    plt.grid(color='#999999',linestyle='dashed')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    
    plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3)
    plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff')
    co0='#555555'
    xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz)
    xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz)
    xs=4.25; ys=dh1; ss='$v_H$=5.0m/s'
    plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz)
    xs=4.25; ys=dh2; ss='$v_H$=2.0m/s'
    plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz)
    co1='black'
    co2='red'
    cont=plt.contour(xx,yy,zh,colors=[co1,co1,co2,co1,co1,co1,co1,co1],levels=[10,15,16.58,20,30,40,60,80])
    cont.clabel(fmt='%1.2f', fontsize=fsz-1)

    plt.subplot(122)
    tstr='Total Cost Contour of Headrace and Penstock (unit: mil.NRP)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Penstock diameter (m)')
    plt.ylabel('Headrace diameter (m)')
    plt.grid(color='#999999',linestyle='dashed')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    
    plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3)
    plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff')
    xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz)
    xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz)
    xs=4.75; ys=dh1; ss='$v_H$=5.0m/s'
    plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz)
    xs=4.75; ys=dh2; ss='$v_H$=2.0m/s'
    plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz)

    cont=plt.contour(xx,yy,zc,colors=[co1],levels=[6000,6500,7000,7500,8000,8200,8500,9000])
    cont.clabel(fmt='%1.0f', fontsize=fsz-1)
    plt.contour(xx,yy,zh,colors=[co2],levels=[16.58])
    xs=5.6; ys=5.65; ss='$h_L$=16.58m'
    plt.text(xs,ys,ss,color='#ff0000',va='bottom',ha='center',rotation=0,fontsize=fsz)
    
    
    fnameF='fig_contour.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def headrace(dd):
    a1=164619 # NRP/m
    b1=450019 # NRP/m
    q=68.64
    al=10149.0
    n=0.016
    v=4*q/np.pi/dd**2
    hl=124.5*n**2*al/dd**(4/3)*v**2/2/9.8
    cc=(a1*dd**2/6.0**2+b1*dd/6.0)*al
    return hl,cc

def penstock(dd):
    a2=143829
    b2=2437167
    cf=177919810
    q=68.64
    al=688.2
    n=0.011
    rho=20.0
    fbr=0.4
    v=4*q/np.pi/dd**2
    h1=124.5*n**2*al/dd**(4/3)
    h2=4*(0.131+0.1632*(dd/rho)**(7/2))
    h3=fbr
    hl=(h1+h2+h3)*v**2/2/9.8+1.991
    cc=(a2*dd**2/4.0**2+b2*dd/4.0)*al+cf
    return hl,cc


def main():
    ds=0.1
    d1=np.arange(4.0,7.0+ds,ds)
    d2=np.arange(3.0,6.0+ds,ds)
    hlt=np.zeros((len(d1),len(d2)),dtype=np.float64)
    cct=np.zeros((len(d1),len(d2)),dtype=np.float64)
    for i,dd1 in enumerate(d1):
        hl1,cc1=headrace(dd1)
        for j,dd2 in enumerate(d2):
            hl2,cc2=penstock(dd2)
            hlt[i,j]=(hl1+hl2)*1.1
            #cct[i,j]=(cc1+cc2)/(450-hlt[i,j])*433.42/1e6
            cct[i,j]=(cc1+cc2)/1e6
            if 16.0<=hlt[i,j]<=16.58: print('{0:6.1f}{1:6.1f}{2:8.3f}{3:8.3f}{4:8.3f}{5:12.3f}'.format(dd1,dd2,hl1*1.1,hl2*1.1,hlt[i,j],cct[i,j]))
    xx, yy = np.meshgrid(d2,d1)
    zh=hlt
    zc=cct
    drawfig(xx,yy,zh,zc)



if __name__ == '__main__':
    main()

以　上

Python 岩盤内埋設式水圧鉄管設計プログラム（旧版）

概要

岩盤内埋設式水圧鉄管の設計プログラムをフルで書いてみたので残しておく。

プログラム前半の非常に長い関数「def xlswrite」は、エクセルに書式指定して結果を書き出すもの。 最終的にワードの報告書に貼り付けなくてはならないので、エクセルにしておくと便利である。このため、面倒ではあったがこの部分も作っておいた。

入力

プログラムでは 'load2.xlsx' を読み込んでいるが、これは出力例「荷重計算」と同じ内容（欄外の変数凡例行は除く）のファイルである。

出力

プログラムでは出力は'penstock2.xlsx' というエクセルファイルに行われ、以下のシートに結果が出力される。

Sheet	Outout
Load	荷重計算（入力ファイルと同じ）
Pin	内圧設計結果
Pex	外圧設計結果

なお設計では、関数「calc1」の中で、水門鉄管技術基準で定める最小板厚とは別に、施工性や耐外厚設計などを考慮し、これより大きな推奨最小板厚（recommended minimum thickness: t0t）というものを導入し、板厚がこれ以下にならないようにしている。

プログラム上での使用鋼材はJIS：SM400, SM490, SM570である。

荷重計算

内圧設計

外圧設計

プログラム

import numpy as np
import pandas as pd
from scipy import optimize
import openpyxl
from openpyxl.styles.borders import Border, Side
from openpyxl.styles import PatternFill
from openpyxl.styles.fonts import Font


def xlline(ws, cmax):
    def wline(ws, row_num, cmax, bl, bc, br):
        col_num=1; ws.cell(row=row_num ,column=col_num).border = bl
        for col_num in range(2,cmax):
            ws.cell(row=row_num ,column=col_num).border = bc
        col_num=cmax; ws.cell(row=row_num ,column=col_num).border = br
    
    # font
    font = Font(name='Courier New',
            size=12,
            bold=False,
            italic=False,
            underline='none',
            strike=False,
            color='000000')
    for row_num in range(1,53):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).font = font
    # cell height and width
    lcol=['A','B','C','D','E','F','G','H','I','J','K','L','M','N']
    for row_num in range(1,53):
        ws.row_dimensions[row_num].height = 20
    for i in range(0,cmax):
        ws.column_dimensions[lcol[i]].width = 12
    ws.column_dimensions[lcol[cmax-1]].width = 30
    
    # set color for cell filling
    fill = PatternFill(patternType='solid', fgColor='ffffaa')
    # write in sheet
    for row_num in range(6,12):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).fill = fill
    for row_num in range(20,26):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).fill = fill
    
    bul = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    buc = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bur = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )
    bcl = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bcc = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bcr = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )
    bll = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    blc = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    blr = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )

    row_num=1; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=2; wline(ws, row_num, cmax, bul, buc, bur)
    for row_num in range(3,29):
        wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=29; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=30; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=31; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=32; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=33; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=34; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=35; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=36; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=37; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=38; wline(ws, row_num, cmax, bll, blc, blr)
    if cmax==11:
        row_num=39; wline(ws, row_num, cmax, bll, blc, blr)
    

def xlswrite(df0,df1,df2):
    # Write results to Excel 
    fnameW='penstock2.xlsx'
    wb = openpyxl.Workbook()
    #
    ws=wb.create_sheet(index=1,title='Load')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='L(m)'
    ws.cell(row=1,column=3).value='Sum(L)'
    ws.cell(row=1,column=4).value='EL(m)'
    ws.cell(row=1,column=5).value='D0(m)'
    ws.cell(row=1,column=6).value='GWL(m)'
    ws.cell(row=1,column=7).value='Hst(m)'
    ws.cell(row=1,column=8).value='Hsg(m)'
    ws.cell(row=1,column=9).value='Hwh(m)'
    ws.cell(row=1,column=10).value='Hin(m)'
    ws.cell(row=1,column=11).value='Hex(m)'
    ws.cell(row=1,column=12).value='Remarks'
    #
    j=1; temp=np.array(df0['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df0['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df0['Sum(L)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df0['EL(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df0['D0(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=np.array(df0['GWL(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df0['Hst(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df0['Hsg(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df0['Hwh(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df0['Hin(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=11; temp=np.array(df0['Hex(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=12; temp=df0['Remarks']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=12; xlline(ws, cmax)
    # legend
    i=38
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'Sum(L): cumulative length'
    i=i+1; ws.cell(row=i,column=1).value = 'EL: elevation of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'GWL: ground water level'
    i=i+1; ws.cell(row=i,column=1).value = 'Hst: hydrostatic head'
    i=i+1; ws.cell(row=i,column=1).value = 'Hsg: pressure head rise by surging'
    i=i+1; ws.cell(row=i,column=1).value = 'Hwh: pressure head rise by water hammering'
    i=i+1; ws.cell(row=i,column=1).value = 'Hin: design internal pressure head'
    i=i+1; ws.cell(row=i,column=1).value = 'Hex: design external pressure head'
    #
    ws=wb.create_sheet(index=1,title='Pin')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='Pi(MPa)'
    ws.cell(row=1,column=3).value='L(m)'
    ws.cell(row=1,column=4).value='D0(mm)'
    ws.cell(row=1,column=5).value='t0(mm)'
    ws.cell(row=1,column=6).value='steel'
    ws.cell(row=1,column=7).value='lam'
    ws.cell(row=1,column=8).value='sig(MPa)'
    ws.cell(row=1,column=9).value='siga(MPa)'
    ws.cell(row=1,column=10).value='Weight(t)'
    ws.cell(row=1,column=11).value='Remarks'
    #
    j=1; temp=np.array(df1['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df1['Pi(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df1['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df1['D0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df1['t0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=df1['steel']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df1['lam'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df1['sig(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df1['siga(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df1['weight(t)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    ws.cell(row=len(temp)+2,column=1).value = 'Sum'
    ws.cell(row=len(temp)+2,column=10).number_format = "0.000" 
    ws.cell(row=len(temp)+2,column=10).value = np.sum(temp)      
    j=11; temp=df0['Remarks']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=11; xlline(ws, cmax)
    # legend
    i=39
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'Pi: design internal pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel'
    i=i+1; ws.cell(row=i,column=1).value = 'lam: internal pressure shared ratio be bedrock'
    i=i+1; ws.cell(row=i,column=1).value = 'sig: stress of pipe shell'
    i=i+1; ws.cell(row=i,column=1).value = 'siga: allowable stress'
    i=i+1; ws.cell(row=i,column=1).value = 'Weight: weight of pipe section'
    #
    ws=wb.create_sheet(index=1,title='Pex')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='Pe(MPa)'
    ws.cell(row=1,column=3).value='L(m)'
    ws.cell(row=1,column=4).value='D0(mm)'
    ws.cell(row=1,column=5).value='t0(mm)'
    ws.cell(row=1,column=6).value='steel'
    ws.cell(row=1,column=7).value='pk_0(MPa)'
    ws.cell(row=1,column=8).value='SF_0'
    ws.cell(row=1,column=9).value='pk_s(MPa)'
    ws.cell(row=1,column=10).value='SF_s'
    ws.cell(row=1,column=11).value='sc_r(MPa)'
    ws.cell(row=1,column=12).value='sc_c(MPa)'
    ws.cell(row=1,column=13).value='SF_c'
    ws.cell(row=1,column=14).value='stiffener'
    #
    j=1; temp=np.array(df2['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df2['Pe(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df2['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df2['D0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df2['t0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=df2['steel']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df2['pk0(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df2['sf0'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df2['pks(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df2['sfs'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=11; temp=np.array(df2['scr(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=12; temp=np.array(df2['scc(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=13; temp=np.array(df2['sfc'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=14; temp=df2['stiffener']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=14; xlline(ws, cmax)
    # legend
    i=38
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'Pe: design external pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel'
    i=i+1; ws.cell(row=i,column=1).value = 'pk_0: critical buckling pressure of pipe section without stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_0: safety factor against external pressure without pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'pk_s: critical buckling pressure of pipe section with stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_s: safety factor against external pressure with stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'sc_r: critical buckling stress of stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'sc_c: compressive stress of stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_c: safety factor against compressive stress in stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'size of stiffener plate: 75mm height x 20mm thickness'
    #    
    wb.save(fnameW)


def stcr(dd0,t0,eps,hr,tr,el,pp,sigF):
    # Calculation of critical stress and compressive stress of stiffener
    def calrg(dd0,t0,eps,hr,tr):
        rm=0.5*(dd0+t0)
        tt=t0-eps
        # change symbols
        b=1.56*np.sqrt(rm*tt)+tr
        d=tt+hr
        s=tt
        t=tr
        # calculation of radius of gyration of area
        e2=(d**2*t+s**2*(b-t))/2/(b*s+t*(d-s))
        e1=d-e2
        ii=1/3*(t*e1**3+b*e2**3-(b-t)*(e2-s)**3)
        aa=s*b+t*(d-s)
        rg=np.sqrt(ii/aa) # radius of gyration of area
        ee=np.abs(e2-s)
        return rg,ee
    def func(x,params):
        # x=sigN
        k0   =params[0]
        rm   =params[1]
        t    =params[2]
        Es   =params[3]
        sigF =params[4]
        rg   =params[5]
        ee   =params[6]
        aa = k0/rm+x/Es
        bb = 1+rm**2/rg**2*x/Es
        cc = 1.68*rm/ee*(sigF-x)/Es*(1-1/4*rm/ee*(sigF-x)/Es)
        f = aa*bb**1.5-cc
        return f
    # main routine
    eps=2.0         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    ns=0.3         # Poisson's ratio of steel plate
    #D0=2100.0  # design internal diameter of penstock
    #t0=30.0    # design plate thickness of penstock
    #sigF=315.0 # yield stress of steel plate (SM490)
    #siga=175.0 # allowable stress of steel plate (SM490)
    t=t0-eps
    rm=(dd0+t0)/2
    r0d=(dd0+2*t0)/2
    k0=0.4e-3*rm
    # calculation of critical stress
    rg,ee=calrg(dd0,t0,eps,hr,tr)
    params=np.array([k0,rm,t,Es,sigF,rg,ee])
    sigN=optimize.brentq(func,0.0,5*sigF,args=(params))
    scr=sigN*(1-r0d/ee*(sigF-sigN)/(1+3/2*np.pi)/Es)
    # calculation of compressive stress
    s0=tr*(t+hr)
    iis=tr/12*(hr+t)**3
    beta=(3*(1-ns**2))**0.25/np.sqrt(rm*t)
    c1=r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) 
    c2=3/(3*(1-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t))
    pd=pp/(tr+1.56*np.sqrt(rm*t))*((tr+1.56*np.sqrt(rm*t))+2*c1/c2)
    sc=pd*r0d*(tr+1.56*np.sqrt(rm*t))/(s0+1.56*t*np.sqrt(rm*t))
    return scr,sc
    

def timo(el,hr,tr,dd0,t0,eps):
    # Calculate criticsl buckling pressure
    # of steel pipe with stiffener by Timoshenko's formula
    # el: interval of stiffener
    # hr: height of stiffener
    # tr: plate thickness of stiffener
    # t: plate thickness
    pi=np.pi
    Es=206000 # elastic modulus of steel
    ns=0.3    # Poisson's ratio of steel
    rm=0.5*(dd0+t0)
    r0d=0.5*(dd0+2*t0)
    t=t0-eps
    s0=tr*(t+hr)
    iis=tr/12*(hr+t)**3
    beta=(3*(1-ns**2))**0.25/np.sqrt(rm*t)
    c1=r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) 
    c2=3/(3*(1-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t))
    lt=2/(tr+1.56*np.sqrt(rm*t))*c1/c2
    lam=1-(1+lt)*(1+tr/(1.56*np.sqrt(rm*t)))/(1+s0/(1.56*t*np.sqrt(rm*t)))
    ell=(el+1.56*np.sqrt(rm*t)*np.arccos(lam))*(1+0.037*np.sqrt(rm*t)/(el+1.56*np.sqrt(rm*t)*np.arccos(lam))*t**3/iis)
    apk=[]
    for n in range(2,30):
        c1=(1-ns**2)/(n**2-1)/(1+n**2*ell**2/pi**2/r0d**2)**2
        c2=t**2/12/r0d**2*((n**2-1)+(2*n**2-1-ns)/(1+n**2*ell**2/pi**2/r0d**2))
        _pk=(c1+c2)*Es*t/(1-ns**2)/r0d
        apk=apk+[_pk]
    pk=min(apk)
    return pk


def ams(dd0,t0,sigF,siga):
    # Calculate criticsl buckling pressure
    # of steel pipe without stiffener by Amstutz's formula
    def func(x,params):
        # x=sigN
        k0   =params[0]
        rm   =params[1]
        t    =params[2]
        Ess  =params[3]
        sigFs=params[4]
        aa = k0/rm+x/Ess
        bb = 1+12*rm**2/t**2*x/Ess
        cc = 3.36*rm/t*(sigFs-x)/Ess*(1-1/2*rm/t*(sigFs-x)/Ess)
        f = aa*bb**1.5-cc
        return f
    def cal_pk(sigN,sigFs,Ess,rm,t):
        aa=rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess)
        pk = sigN/aa
        return pk
    def cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag):
        aa = alphas*deltaT+betag*siga*eta/Es
        k0 = aa*r0d/(1+betag)
        return k0    
    # main routine
    eps=2.0         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    pos=0.3         # Poisson's ratio of steel plate
    eta=0.85         # welding efficiency
    alphas=1.2e-5   # thermal expansion coefficient
    deltaT=20.0     # temperature decrease
    betag=1.0       # plastic deformation coefficient of bedrock
    #D0=2100.0  # design internal diameter of penstock
    #t0=30.0    # design plate thickness of penstock
    #sigF=315.0 # yield stress of steel plate (SM490)
    #siga=175.0 # allowable stress of steel plate (SM490)
    Ess=Es/(1-pos**2)
    mu=1.5-0.5*1/(1+0.002*Es/sigF)**2
    sigFs=mu*sigF/(1-pos+pos**2)**0.5
    t=t0-eps
    rm=(dd0+t0)/2
    r0d=(dd0+2*t0)/2
    k0=cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag)
    k0=0.4e-3*rm
    params=np.array([k0,rm,t,Ess,sigFs])
    sigN=optimize.brentq(func,0.0,sigFs,args=(params))
    pk=cal_pk(sigN,sigFs,Ess,rm,t)
    return pk


def datainp():
    df0 = pd.read_excel('load2.xlsx')
    return df0


def calc1(num,pp,dd0,sa,eps):
    # Design for internal pressure
    # Calculate required plate thickness
    Es=206000.0
    alpha=1.2e-5
    deltaT=20.0
    Ec=20600.0
    Bc=0.0
    Eg=2000.0
    Bg=1.0
    mg=4.0
    dr=5.2
    
    t0t=25 # recommended minimum thickness
    t0min=np.ceil((dd0+800)/400)
    dd=dd0-eps
    lam=0
    tt=pp*dd/2/sa
    t0=np.ceil(tt+eps)
    if t0<t0min: t0=t0min
    if t0<t0t: t0=t0t
    # Internal pressure shared design by bedrock
    if 5<=num<=10 or 20<=num<=25:
        c1=sa-Es*alpha*deltaT
        c2=(1+Bc)*Es/Ec*2/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2/dd
        tt=pp*dd/2/sa-c1/c2/sa
        lam1=1-Es/pp*alpha*deltaT*2*tt/dd
        lam2=1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd
        lam=lam1/lam2
        t0=np.ceil(tt+eps)
        if t0<t0min: t0=t0min
        if t0<t0t: t0=t0t
        tt=t0-eps
        lam1=1-Es/pp*alpha*deltaT*2*tt/dd
        lam2=1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd
        lam=lam1/lam2    
    sig=pp*dd/2/(t0-eps)*(1-lam)
    return t0,sig,lam


def main():    
    df0=datainp()
    no=np.array(df0['No'],dtype=np.int)
    al=np.array(df0['L(m)'],dtype=np.float64)        # length (m)
    d0=np.array(df0['D0(m)'],dtype=np.float64)*1000  # internal diameter (mm)
    pi=np.array(df0['Hin(m)'],dtype=np.float64)*0.01 # internal pressure (MPa)
    pe=np.array(df0['Hex(m)'],dtype=np.float64)*0.01 # external pressure (MPa)
    eps=2.0
    ef=0.85

    a_t0=np.zeros(len(no),dtype=np.float64)
    a_sig=np.zeros(len(no),dtype=np.float64)
    a_lam=np.zeros(len(no),dtype=np.float64)
    a_sa=np.zeros(len(no),dtype=np.float64)
    a_sname=[]
    a_pk0=np.zeros(len(no),dtype=np.float64)
    a_sf0=np.zeros(len(no),dtype=np.float64)
    a_pks=np.zeros(len(no),dtype=np.float64)
    a_sfs=np.zeros(len(no),dtype=np.float64)
    a_scr=np.zeros(len(no),dtype=np.float64)
    a_scc=np.zeros(len(no),dtype=np.float64)
    a_sfc=np.zeros(len(no),dtype=np.float64)
    a_stiff=[]
    
    for num,pp,dd0,ppe in zip(no,pi,d0,pe):
        sname='SM400'
        sa=130*ef
        t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM400' and 40<t0:
            sname='SM490'
            sa=175*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM490' and 40<t0:
            sname='SM570'
            sa=240*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM570' and 40<t0:
            sname='SM570'
            sa=235*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        a_t0[num-1]=t0
        a_sig[num-1]=sig/ef
        a_lam[num-1]=lam
        a_sa[num-1]=sa/ef
        a_sname=a_sname+[sname]
        if sname=='SM400': sigF=235
        if sname=='SM490': sigF=315
        if sname=='SM570' and t0<=40: sigF=450
        if sname=='SM400' and 40<t0: sigF=430
        # without stiffener
        pk0=ams(dd0,t0,sigF,sa/ef)
        sf0=pk0/ppe
        if sf0<1.5:
            # stiffener
            hr=75; tr=20
            # stiffener interval=2000mm
            el=3000.0
            pk1=timo(el,hr,tr,dd0,t0,eps)
            sf1=pk1/ppe
            scr1,sc1=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc1=scr1/sc1
            # stiffener interval=1500mm
            el=1500.0
            pk2=timo(el,hr,tr,dd0,t0,eps)
            sf2=pk2/ppe
            scr2,sc2=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc2=scr2/sc2
            # stiffener interval=1000mm
            el=1000.0
            pk3=timo(el,hr,tr,dd0,t0,eps)
            sf3=pk3/ppe
            scr3,sc3=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc3=scr3/sc3
        # store results
        a_pk0[num-1]=pk0
        a_sf0[num-1]=sf0        
        if 1.5<=sf0:
            a_pks[num-1]=0
            a_sfs[num-1]=0
            a_stiff=a_stiff+['no-need']
            a_scr[num-1]=0
            a_scc[num-1]=0
            a_sfc[num-1]=0
        elif sf0<1.5 and 1.5<=sf1:
            a_pks[num-1]=pk1
            a_sfs[num-1]=sf1
            a_stiff=a_stiff+['interval=3000mm']
            a_scr[num-1]=scr1
            a_scc[num-1]=sc1
            a_sfc[num-1]=sfc1
        elif sf1<1.5 and 1.5<=sf2:
            a_pks[num-1]=pk2
            a_sfs[num-1]=sf2
            a_stiff=a_stiff+['interval=1500mm']
            a_scr[num-1]=scr2
            a_scc[num-1]=sc2
            a_sfc[num-1]=sfc2
        elif sf2<1.5 and 1.5<=sf3:
            a_pks[num-1]=pk3
            a_sfs[num-1]=sf3
            a_stiff=a_stiff+['interval=1000mm']
            a_scr[num-1]=scr3
            a_scc[num-1]=sc3
            a_sfc[num-1]=sfc3
    
    ww=0.25*np.pi*((d0+2*a_t0)**2-d0**2)/1000/1000*al*7.85 # weight
    print('Weight:{0:10.3f} ton'.format(np.sum(ww)))
    print('Average thickness: {0:6.1f} mm'.format(np.sum(a_t0*al)/np.sum(al)))
    
    df1=pd.DataFrame({
        'No': no,
        'Pi(MPa)': pi,
        'L(m)' : al,
        'D0(mm)': d0,
        't0(mm)': a_t0,
        'steel': a_sname,
        'lam': a_lam,
        'sig(MPa)': a_sig,
        'siga(MPa)': a_sa,
        'weight(t)': ww,
        'Remarks': df0['Remarks']
    })
    df2=pd.DataFrame({
        'No': no,
        'Pe(MPa)': pe,
        'L(m)' : al,
        'D0(mm)': d0,
        't0(mm)': a_t0,
        'steel': a_sname,
        'pk0(MPa)': a_pk0,
        'sf0': a_sf0,
        'pks(MPa)': a_pks,
        'sfs': a_sfs,
        'scr(MPa)': a_scr,
        'scc(MPa)': a_scc,
        'sfc': a_sfc,
        'stiffener': a_stiff
    })
    with pd.ExcelWriter('penstock.xlsx') as writer:
        df0.to_excel(writer,sheet_name='Load',index=False)
        df1.to_excel(writer,sheet_name='Internal Pressure',index=False)
        df2.to_excel(writer,sheet_name='External Pressure',index=False)

        
    xlswrite(df0,df1,df2)
        
        
if __name__ == '__main__':
    main()

以　上

Python jpg画像を縮小保存しhtmlで表示

タイトル通り、jpg画像を縮小保存しhtmlで表示する。

import glob
import os
from PIL import Image
import os

def select_pic():
    path = '*.jpg'
    # 現フォルダ内にあって拡張子がJPGのファイル名を取得
    file_list = glob.glob(path)
    print(file_list)
    return file_list

def process(file_list):
    print(file_list)
    #for frameR in file_list:
    for i in range(len(file_list)):
        fnameR=file_list[i]
        print(fnameR)
        # 元画像の読み込み
        originalImg = Image.open(fnameR)
        # 元画像のサイズを取得
        width1, height1 = originalImg.size
        # サムネイル作成
        width2=600
        height2=width2*height1/width1
        originalImg.thumbnail((int(width2), int(height2)), Image.ANTIALIAS)
        name=fnameR.replace('.jpg','')
        fnameW='600_'+name+'.jpg'
        # 保存
        originalImg.save(fnameW)


def make_html():
    filenames = os.listdir('./')
    imgl=[]
    ww=[]
    hh=[]
    for fname in sorted(filenames):
        if fname[0:4]=='600_':
            print(fname)
            im=Image.open(fname)
            w=im.size[0]
            h=im.size[1]
            imgl=imgl+[fname]
            ww=ww+[w]
            hh=hh+[h]
    f=open('maggie.html','w')
    print('<html>',file=f)
    print('<body>',file=f)
    print('<table>',file=f)
    n=len(imgl)
    m=int(n/5)+1
    k=-1
    for i in range(0,m):
        print('<tr>',file=f)
        for j in range(0,5):
            k=k+1
            if k<=n-1:
                pic=imgl[k]
                w1=200
                h1=int(hh[k]/ww[k]*200)
                print('<td align="center"><img src="'+pic+'" alt="pic" width="'+str(w1)+'", height="'+str(h1)+'"><br><a href="'+pic+'">'+pic+'<a></td>',file=f)
            else:
                print('<td></td>',file=f)
        print('</tr>',file=f)
    print('</table>',file=f)
    print('</body>',file=f)
    print('</html>',file=f)
    f.close()

def main():
    file_list=select_pic()
    process(file_list)
    make_html()


if __name__ == '__main__':
    main()

以　上

雑記 また機械学習の本を買ってしまった

記事の最後に行く

昨日（2019.08.11）、また機械学習の本を買ってしまった。

• オライリー・ジャパン　PythonによりAIプログラミング入門　ディープラーニングを始める前に身に着けておくべき15の基礎技術

あと、久しぶりにpipで一括アップデートしたら、matplotlibがなんとなくおかしい。

前にjupyterがおかしかったときは、tornadoの問題だった（以下参照）

• https://qiita.com/Kit-Ok/items/077a4033057ba4aaf990

今回はmatplotlibの問題のようです。

• https://qiita.com/yoshi65/items/90532732bf9d3875bec7

Thank you.

記事の先頭に行く

雑記 機械学習の勉強を始めてみる

記事の最後に行く

いよいよ７月。 昨年（2018年10月10日）、マレーシアでの業務が完了し、帰国した。 それ以降、以下のような感じで過ごしてきた。

• 2018年10月 ひま
• 2018年11月 ひま
• 2018年12月 忙しい
• 2019年01月 忙しい
• 2019年02月 忙しい（ベトナム出張２週間）
• 2019年03月 忙しい
• 2019年04月 上旬は忙しい（マレーシア出張２週間）、中旬以降ひま
• 2019年05月 ひま
• 2019年06月 ひま

昨年の10月11月は急速の意味もありひまでよかったのだが、今年になって4月中旬以降の2ヶ月半はなんだか時間を無駄に過ごしてきた感じ。

そこで昨日久しぶりに本を買った。

• オライリー・ジャパン「Pythonではじめる機械学習」

scikit-learn を用いて機械学習の基礎を勉強するもの。

機械学習もちらほらネットを見ながらかじっては見ていたが、やはり本に沿って一式流してみないと体系的な知識がつかないなと感じていたところ。

今の予定では８月からまた忙しくなる「予定」ではあるが、第１章を読んだところ、比較的わかりやすく書かれているようなので、なんとか最後まで読破しようと思っている。

ついでに英語の勉強も少しやってみようかな。

Thank you.

記事の先頭に行く

設計 Pythonによる水文量頻度解析用自作関数とその利用

記事の最後に行く

自作関数

以下の確率分布に対するパラメータ推定と予測値計算の部分を自作関数としてプログラム本体から分離している。

含まれる関数（確率分布）

関数名	説明
ln3	３変数対数正規分布（LN3）
lp3	対数ピアソンIII型分布（LP3）
gev	一般化極値分布（GEV）
gum	グンベル分布（GUM）
sqe	平方根指数型最大値分布（SQRT-ET）

関数への入力値（引数）はいずれについても以下の３つ。

入力値（引数）

変数名	説明
xx	昇順に並び替えられた水文量観測値
pp	各観測値に相当する非超過確率（プロッティング・ポジション公式より算出）<.td>
pyear	予測値を求めたい確率年（昇順）

関数からの戻り値はいずれについても以下の４つ。

戻り値

変数名	説明
res	計算された確率分布関数のパラメータ（numpy.array）
ye	入力した非超過確率（pp）に対応する予測値（numpy.array）
yp	入力した確率年に対する予測値（numpy.array）
ss	確率分布名、算定されたパラメータなどの情報を記載した文字列

#im_hydrology.py
import numpy as np
from scipy.stats import norm  # normal distribution for LN3
from scipy.stats import gamma # Gamma distribution for LP3
import scipy.special as ssp   # Gamma function for GEV

def ln3(xx,pp,pyear):
    def func_ln3(xx):
        n=len(xx)
        mx=np.sum(xx)/n
        s1=np.sqrt(np.sum((xx-mx)**2)/n)
        cx=np.sum(((xx-mx)/s1)**3)/n
        sx=np.sqrt(n/(n-1))*s1
        #rx=np.sqrt(n*(n-1))/(n-2)*cx
        a=1.01+7.01/n+14.66/n**2
        b=1.69/n+74.66/n**2
        rx=cx*(a+b*cx**3)
        beta=1+rx**2/2
        x=(beta+np.sqrt(beta**2-1))**(1/3)+(beta-np.sqrt(beta**2-1))**(1/3)-1
        sd=np.sqrt(np.log(x))
        mu=np.log(sx/np.sqrt(x*(x-1)))
        aa=mx-np.exp(mu)*np.exp(sd**2/2)
        return aa,mu,sd
    def est_ln3(aa,mu,sd,pp):
        zz=norm.ppf(pp, loc=0, scale=1)
        return aa+np.exp(mu+sd*zz)
    # parameter calculation
    aa,mu,sd=func_ln3(xx)
    # estimation
    ye=est_ln3(aa,mu,sd,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_ln3(aa,mu,sd,prob)
    res=np.array([aa,mu,sd,rr])
    ss='LN3\na={0:8.3f}\nm={1:8.3f}\ns={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def lp3(xx,pp,pyear):
    def func_lp3(xx):
        n=len(xx)
        yy=np.log(xx)
        my=np.sum(yy)/n
        sy=np.sqrt(np.sum((yy-my)**2)/n)
        cy=np.sum(((yy-my)/sy)**3)/n
        sy=np.sqrt(n/(n-1))*sy
        a=1+6.51/n+20.2/n**2
        b=1.48/n+6.77/n**2
        ry=cy*(a+b*cy**2)
        bb=4/ry**2
        aa=sy/np.sqrt(bb)
        if ry<0: aa=-aa
        cc=my-aa*bb
        return bb,cc,aa
    def est_lp3(bb,cc,aa,pp):
        ww=gamma.ppf(pp, bb, loc=0, scale=1)
        return np.exp(cc+aa*ww)
    # parameter calculation
    bb,cc,aa=func_lp3(xx)
    # estimation
    ye=est_lp3(bb,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_lp3(bb,cc,aa,prob)
    res=np.array([bb,cc,aa,rr])
    ss='LP3\nb={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gev(xx,pp,pyear):
    def func_gev(xx):
        n=len(xx)
        jj=np.arange(n)+1
        b0=np.sum(xx)/n
        b1=np.sum((jj-1)*xx)/n/(n-1)
        b2=np.sum((jj-1)*(jj-2)*xx)/n/(n-1)/(n-2)
        lam1=b0
        lam2=2*b1-b0
        lam3=6*b2-6*b1+b0
        d=2*lam2/(lam3+3*lam2)-np.log(2)/np.log(3)
        k=7.8590*d+2.9554*d**2
        a=k*lam2/(1-2**(-k))/ssp.gamma(1+k)
        c=lam1-a/k*(1-ssp.gamma(1+k))
        return k,c,a
    def est_gev(kk,cc,aa,pp):
        return cc+aa/kk*(1-(-np.log(pp))**kk)
    # parameter calculation
    kk,cc,aa=func_gev(xx)
    # estimation
    ye=est_gev(kk,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gev(kk,cc,aa,prob)
    res=np.array([kk,cc,aa,rr])
    ss='GEV\nk={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gum(xx,pp,pyear):
    def func_gum(xx):
        n=len(xx)
        b0=np.sum(xx)/n
        j=np.arange(0,n)
        b1=np.sum(j*xx)/n/(n-1)
        lam1=b0
        lam2=2*b1-b0
        aa=lam2/np.log(2)
        cc=lam1-0.5772*aa
        return aa,cc
    def est_gum(aa,cc,pp):
        return cc-aa*np.log(-np.log(pp))
    aa,cc=func_gum(xx)
    # estimation
    ye=est_gum(aa,cc,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gum(aa,cc,prob)
    res=np.array([aa,cc,rr])
    ss='GUM\na={0:8.3f}\nc={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss


def sqe(xx,pp,pyear):
    def est_sqrtet(aa,bb,pp,x0):
        # estimation of quantile by Newton-Raphson Method
        def gx(x,aa,p):
            ggx=np.log(1+x)-x-np.log(-1/aa*np.log(p))
            gdx=1/(1+x)-1
            return ggx,gdx
        imax=100
        eps=1.0e-10
        tp=np.zeros(len(pp),dtype=np.float64)
        x1=x0
        for i,p in enumerate(pp):
            for j in range(imax):
                ggx,gdx=gx(x1,aa,p)
                x2=x1-ggx/gdx
                if np.abs(x2-x1)<eps: break
                x1=x2
            tp[i]=x2
        return tp**2/bb

    def bisection(xx):
        # estimation of coefficient bb by Bisection Method
        def func(bb,xx):
            # function for Bisection Method
            nn=len(xx)
            a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
            a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx))))
            return a1-a2
        imax=100    # maximum number of iterations
        err=1.0e-10 # allowable error of the zero value
        x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2
        x2=x1+0.5
        for i in range(0,imax):
            xi=0.5*(x1+x2)
            f1=func(x1,xx)
            f2=func(x2,xx)
            fi=func(xi,xx)
            if 0<f1*f2:
                x1=x2
                x2=x1+0.5
            elif abs(x2-x1) < err:
                break # converged
            else:
                if f1*fi < 0.0: x2=xi
                if fi*f2 < 0.0: x1=xi
        return xi

    bb=bisection(xx)
    aa=(np.sum(np.sqrt(bb*xx))-2*len(xx))/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
    # estimation
    x0=np.sqrt(bb*xx[-1])
    ye=est_sqrtet(aa,bb,pp,x0)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_sqrtet(aa,bb,prob,x0)
    res=np.array([aa,bb,rr])
    ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss

クオンタイルプロット

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology


fsz=12
xmin=0; xmax=5000; dx=1000
ymin=0; ymax=5000; dy=1000


def qqplot(fnameR,xx,ye,ss,n):
    tstr=fnameR
    nnn=230+n
    plt.subplot(nnn)
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Observed (m$^3$/s)')
    plt.ylabel('Predicted (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.plot(xx,ye,'o',clip_on=False)
    plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1)
    xs=xmin+0.05*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished')


def calc(fnameR,alpha,pyear):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    res1,ye1,yp1,ss1=im_hydrology.ln3(xx,pp,pyear)    
    res2,ye2,yp2,ss2=im_hydrology.lp3(xx,pp,pyear)    
    res3,ye3,yp3,ss3=im_hydrology.gev(xx,pp,pyear)
    res4,ye4,yp4,ss4=im_hydrology.gum(xx,pp,pyear)
    res5,ye5,yp5,ss5=im_hydrology.sqe(xx,pp,pyear)
    
    plt.figure(figsize=(12,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    qqplot(fnameR,xx,ye1,ss1,1)
    qqplot(fnameR,xx,ye2,ss2,2)
    qqplot(fnameR,xx,ye3,ss3,3)
    qqplot(fnameR,xx,ye4,ss4,4)
    qqplot(fnameR,xx,ye5,ss5,5)
    plt.tight_layout()
    s=fnameR.replace('inp_','')
    s=s.replace('.txt','')
    fnameF='fig_qq_'+s+'.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()    

    s0='{0:8s}'.format('Year')
    s1='{0:8s}'.format('LN3')
    s2='{0:8s}'.format('LP3')
    s3='{0:8s}'.format('GEV')
    s4='{0:8s}'.format('GUM')
    s5='{0:8s}'.format('SQRT-ET')
    for i in range(len(pyear)):
        s0=s0+',{0:6.0f}'.format(pyear[i])
        s1=s1+',{0:6.0f}'.format(yp1[i])
        s2=s2+',{0:6.0f}'.format(yp2[i])
        s3=s3+',{0:6.0f}'.format(yp3[i])
        s4=s4+',{0:6.0f}'.format(yp4[i])
        s5=s5+',{0:6.0f}'.format(yp5[i])
    print(s0)
    print(s1)
    print(s2)
    print(s3)
    print(s4)
    print(s5)


def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        calc(fnameR,alpha,pyear)


#==============
# Execution
#==============
if __name__ == '__main__': main()

確率年と予測値

# Jackknife method
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology

# Global variables
fsz=14
xmin=1
xmax=1000
dx=1
ymin=0
ymax=10000
dy=1000


def drawfig0(pyear,qqq,smt,ss):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qqq[0,:],'o--',ms=10,color=col,markerfacecolor=col,label=smt[0],lw=1,clip_on=False)
    plt.plot(pyear,qqq[1,:],'s--',ms=10,color=col,markerfacecolor=col,label=smt[1],lw=1,clip_on=False)
    plt.plot(pyear,qqq[2,:],'^--',ms=10,color=col,markerfacecolor=col,label=smt[2],lw=1,clip_on=False)
    plt.plot(pyear,qqq[3,:],'o-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[3],lw=1,clip_on=False)
    plt.plot(pyear,qqq[4,:],'s-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[4],lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    fnameF='fig_comp.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def drawfig(pyear,qb,sd,ss,kkk):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    if kkk==1: tstr=tstr+' (LN3)'
    if kkk==2: tstr=tstr+' (LP3)'
    if kkk==3: tstr=tstr+' (GEV)'
    if kkk==4: tstr=tstr+' (GUM)'
    if kkk==5: tstr=tstr+' (SQRT-ET)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qb   ,'o-'  ,ms=10,color=col,markerfacecolor=col,label='Average',lw=1,clip_on=False)
    plt.plot(pyear,qb+sd,'^--' ,ms=10,color=col,markerfacecolor=col,label='Ave.+ SD',lw=1,clip_on=False)
    plt.plot(pyear,qb-sd,'s--' ,ms=10,color=col,markerfacecolor=col,label='Ave.- SD',lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    if kkk==1: fnameF='fig_rp_ln3.jpg'
    if kkk==2: fnameF='fig_rp_lp3.jpg'
    if kkk==3: fnameF='fig_rp_gev.jpg'
    if kkk==4: fnameF='fig_rp_gum.jpg'
    if kkk==5: fnameF='fig_rp_sqe.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def calc(fnameR,alpha,pyear,kkk):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    if kkk==1: res,ye0,yp0,ss=im_hydrology.ln3(xx,pp,pyear)    
    if kkk==2: res,ye0,yp0,ss=im_hydrology.lp3(xx,pp,pyear)    
    if kkk==3: res,ye0,yp0,ss=im_hydrology.gev(xx,pp,pyear)    
    if kkk==4: res,ye0,yp0,ss=im_hydrology.gum(xx,pp,pyear)
    if kkk==5: res,ye0,yp0,ss=im_hydrology.sqe(xx,pp,pyear)
    print(ss)
    
    qq=np.zeros((n,len(pyear)),dtype=np.float64)
    qm=np.zeros(len(pyear),dtype=np.float64)
    qb=np.zeros(len(pyear),dtype=np.float64)
    sd=np.zeros(len(pyear),dtype=np.float64)
    for i in range(0,n):
        xi=np.delete(xx,i)
        m=len(xi)
        ii=np.arange(m)+1
        pp=(ii-alpha)/(m+1-2*alpha)
        if kkk==1: res,ye,yp,ss=im_hydrology.ln3(xi,pp,pyear)
        if kkk==2: res,ye,yp,ss=im_hydrology.lp3(xi,pp,pyear)
        if kkk==3: res,ye,yp,ss=im_hydrology.gev(xi,pp,pyear)
        if kkk==4: res,ye,yp,ss=im_hydrology.gum(xi,pp,pyear)
        if kkk==5: res,ye,yp,ss=im_hydrology.sqe(xi,pp,pyear)
        qq[i,:]=yp
    for j in range(0,len(pyear)):
        qm[j]=np.average(qq[:,j])
        qb[j]=n*yp0[j]-(n-1)*qm[j]
        sd[j]=np.sqrt((n-1)/n*np.sum((qq[:,j]-qm[j])**2))
    return qb,sd
    

def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    qqq=np.zeros((5,len(pyear)),dtype=np.float64)
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        for kkk in [1,2,3,4,5]:
            qb,sd=calc(fnameR,alpha,pyear,kkk)
            s0='{0:4s}'.format('Year')
            s1='{0:4s}'.format('Ave.')
            s2='{0:4s}'.format('SD')
            for j in range(len(pyear)):
                s0=s0+'{0:6.0f}'.format(pyear[j])
                s1=s1+'{0:6.0f}'.format(qb[j])
                s2=s2+'{0:6.0f}'.format(sd[j])
            ss='Retuen Period and Flood Discharge (m$^3$/s)'
            ss=ss+'\n'+s0+'\n'+s1+'\n'+s2
            ss=ss+'\n* Ave.: average, SD: standard deviation'
            print(ss)
            drawfig(pyear,qb,sd,ss,kkk)
            qqq[kkk-1,:]=qb[:]
        smt=['LN3','LP3','GEV','GUM','SQRT-ET']
        s0='{0:8s}'.format('Year')
        s1='{0:8s}'.format(smt[0])
        s2='{0:8s}'.format(smt[1])
        s3='{0:8s}'.format(smt[2])
        s4='{0:8s}'.format(smt[3])
        s5='{0:8s}'.format(smt[4])
        for j in range(len(pyear)):
            s0=s0+'{0:6.0f}'.format(pyear[j])
            s1=s1+'{0:6.0f}'.format(qqq[0,j])
            s2=s2+'{0:6.0f}'.format(qqq[1,j])
            s3=s3+'{0:6.0f}'.format(qqq[2,j])
            s4=s4+'{0:6.0f}'.format(qqq[3,j])
            s5=s5+'{0:6.0f}'.format(qqq[4,j])
        ss='Retuen Period and Flood Discharge (m$^3$/s)'
        ss=ss+'\n'+s0+'\n'+s1+'\n'+s2+'\n'+s3+'\n'+s4+'\n'+s5
        print(ss)
        drawfig0(pyear,qqq,smt,ss)


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く

FEM 格子桁構造解析プログラム

記事の最後に行く

はじめに

格子桁構造解析プログラムを紹介します。 このプログラムは、平面骨組構造解析プログラムの軸力・軸方向変位をねじりモーメント・ねじり角に書き換えたものです。 構造は平面形状とし、x - y 平面上で定義します。

格子桁解析は、立体骨組構造解析プログラムで行うことができるのですが、床あるいは壁のような単純な板状構造物の解析には、入力データの単純さを考えると格子桁解析も有用であることから、作ったものです。

なお、このプログラムでは温度変化による断面力は考慮できません。 それは、考慮できる断面力は、ねじりモーメント、曲げモーメント、せん断力であり、軸力項が無いためです。

（注意事項）

充実矩形断面のねじり定数 は、断面の短辺を 、長辺を として、以下の式で計算できます。

ここで、コンクリート床のような版構造を格子桁でモデル化する場合、ねじり定数 はゼロとみなせるくらいに小さな値を入力します。 これは、ある要素のねじり角は、それに交差する要素の曲げ角に影響するため、曲げモーメントがねじりの影響で不連続になってしまうためです。

本来の格子桁をモデル化する場合は、もちろん、計算により求めた正規のねじり定数を入力する必要があります。

Grid Girder Analysis

Outline of Grid Girder Analysis Program

• This is a program for Grid Girder Analysis. Elastic and small displacement problems can be treated. This program has been made using 'Python 3.'
• Beam element with 2 nodes is used. 1 node has 3 degrees of freedom which are torsional rotation, bending rotation and deflection.
• Simultaneous linear equations are solved using numpy.linalg.solve.
• Input / Output file name can be defined arbitrarily as text files with the separator of blank, and those are inputted from command line of a terminal.

Workable condition

Item	Description
Element	Beam element which has 3 degrees of freedom such as torsional rotation, bending rotation and deflection.
Load	Specify the node number and load values on the node.
Displacement	Specify the node number and displacements on the node. Any values can be applied including zero

Global Coordinate System	Local Coordinate Syatem

The Stiffness matrix of an element is shown below:

A shear modulus is calculated from an elastic modulus of and Poisson's ration using following equation in the program.

Since the coordinate transformation is carried out on the x-y plane only, the transfoemation matrix is the same type as it for 2D frame analysis.

Format of input data file

npoin nele nsec npfix nlod  # Basic values for analysis
E po AA AI AJ gamma         # material properties
    ....(1 to nsec)....     #
node-1,node-2,matno         # Element connectivity, material set number, uniformly distributed load per unit length
    ....(1 to nele)....     #
x,y                         # Node coordinates, temperature change of node
    ....(1 to npoin)....    #
node nokx noky nokz rdisx rdisy rdisz # node number, restrict conditions and forced displacement
    ....(1 to npfix)....    #
node,fx,fy,fz               # node number, Load values of torsional moment, bending moment, vertical force
    ....(1 to nlod)....     # (Omit data input if nlod=0)

npoin : Number of nodes                   E     : Elastic modulus of element
nele  : Number of elements                po    : Poison's ratio
nsec  : Number of material sets           AA    : Section area of element
npfix : Number of restricted nodes        AI    : Moment of Inertia of element
nlod  : Number of loaded nodes            AJ    : Tosion constant
                                          gamma : unit weight of material
                                          matno : Material set number

• The structure shall be defined on the x-y plane in the global coordinate system.
• x-direction in the local coordinate system means the member axis.
• The displacement in local x-direction means the rotation around the x-axis, which is same as torsional rotation.
• The displacement in local y-direction means the rotation around the y-axis, which is same as bending rotation.
• The displacement in local z-direction means the deflection of the beam.
• Restricted node means the node which has known (given) displacement. As a known (given) value of nodal displacement, any value can be given including zero for a restricted node.
• Since stress resultants of element are defined as equivalent nodal forces in local coordinate system, it is necessary to note that signs are different from it on general structural mechanics. Positive directions are always right-direction, upward-direction for each node in local coordinate system.

Format of output file

npoin  nele  nsec npfix  nlod
    (Each value for above item)
  sec   E   po   A   I   J   gamma
    (material caracteristics)
    .....(1 to nsec).....
 node   x   y   fx   fy   fz   kox   koy   koz
    node: node number
    x   : x-coordinate
    y   : y-coordinate
    fx  : torsional load
    fy  : bending moment load
    fz  : vertical force load
    kox : restriction around x-axis
    koy : restriction around y-axis
    koz : restriction in z-direction
    .....(1 to npoin).....
 node   kox   koy   koz   rdis_x   rdis_y   rdis_z
    node   : node number
    kox    : restriction around x-axis
    koy    : restriction around y-axis
    koz    : restriction in z-direction
    rdis_x : given rotation around x-axis
    rdis_y : given rotation around y-axis
    rdis_z : given displacement in z-direction
    .....(1 to npfix).....
 elem   i   j   sec
    elem : element number
    i    : node-1
    j    : node-2
    sec  : material number
    .....(1 to nele).....
 node   rot-x   rot-y   dis-z
    node  : node number
    rot-x : rotation arouns x-axis
    rot-y : rotation around y-axis
    dis-z : displacement in z-direction
    .....(1 to npoin).....
 elem   T_i   M_i   Q_i   T_j   M_j   Q_j
    elem : element number
    T_i  : torsional moment at node-i
    M_i  : bending moment at node-i
    Q_i  : shearing force at node-i
    T_j  : torsional moment at node-j
    M_j  : bending moment at node-j
    Q_j  : shearing force at node-j
    .....(1 to nele).....
n=693  time=0.078 sec   # reference information

Grid Girder Analysis Program

###################################
# Grid Girder Analysis
###################################
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_grd(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    # array declaration
    x     =np.zeros([2,npoin],dtype=np.float64)         # Coordinates of nodes
    ae    =np.zeros([6,nsec],dtype=np.float64)          # Section characteristics
    node  =np.zeros([nod+1,nele],dtype=np.int)          # Node-element relationship
    fp    =np.zeros(nfree*npoin,dtype=np.float64)       # External force vector
    mpfix =np.zeros([nfree,npoin],dtype=np.int)         # Ristrict conditions
    rdis  =np.zeros([nfree,npoin],dtype=np.float64)     # Ristricted displacement
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) # E  : Elastic modulus
        ae[1,i]=float(text[1]) # po : Poisson's ratio
        ae[2,i]=float(text[2]) # AA : Section area
        ae[3,i]=float(text[3]) # AI : Moment of inertia
        ae[4,i]=float(text[4]) # AJ : Tortional constant
        ae[5,i]=float(text[5]) # gamma : Unit weight of beam material
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]=float(text[0])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])             #fixed node
        mpfix[0,lp-1]=int(text[1])  #index of restriction for tortional rotation
        mpfix[1,lp-1]=int(text[2])  #index of restriction for bending rotation
        mpfix[2,lp-1]=int(text[3])  #index of restriction for vertical displacement
        rdis[0,lp-1]=float(text[4]) #given tortional rotation
        rdis[1,lp-1]=float(text[5]) #given bending rotation
        rdis[2,lp-1]=float(text[6]) #given vertical displacement
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[3*lp-3]=float(text[1]) #tortional moment load
            fp[3*lp-2]=float(text[2]) #bending moment load
            fp[3*lp-1]=float(text[3]) #vertical load
    f.close()
    return npoin, nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp


def prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('sec','E','po','A','I','J','gamma'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
    .format('node','x','y','fx','fy','fz','kox','koy','koz'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d} {8:5d}'
        .format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>15s} {5:>15s} {6:>15s}'
    .format('node','kox','koy','koz','rdis_x','rdis_y','rdis_z'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]+mpfix[2,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:15.7e} {5:15.7e} {6:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],rdis[0,i],rdis[1,i],rdis[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

    
def prout_grd(fnameW,npoin,nele,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','rot-x','rot-y','dis-z'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','T_i','M_i','Q_i','T_j','M_j','Q_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
    fout.close()


def sm_grd(GJ,EI,x1,y1,x2,y2):
    ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
    tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
    xx=x2-x1
    yy=y2-y1
    el=np.sqrt(xx**2+yy**2)
    ek[0,0]= GJ/el; ek[0,1]=        0.0; ek[0,2]=         0.0; ek[0,3]=-GJ/el; ek[0,4]=        0.0; ek[0,5]=         0.0
    ek[1,0]=   0.0; ek[1,1]= 4*EI/el   ; ek[1,2]= -6*EI/el**2; ek[1,3]=   0.0; ek[1,4]= 2*EI/el   ; ek[1,5]=  6*EI/el**2
    ek[2,0]=   0.0; ek[2,1]=-6*EI/el**2; ek[2,2]= 12*EI/el**3; ek[2,3]=   0.0; ek[2,4]=-6*EI/el**2; ek[2,5]=-12*EI/el**3
    ek[3,0]=-GJ/el; ek[3,1]=        0.0; ek[3,2]=         0.0; ek[3,3]= GJ/el; ek[3,4]=        0.0; ek[3,5]=         0.0
    ek[4,0]=   0.0; ek[4,1]= 2*EI/el   ; ek[4,2]= -6*EI/el**2; ek[4,3]=   0.0; ek[4,4]= 4*EI/el   ; ek[4,5]=  6*EI/el**2
    ek[5,0]=   0.0; ek[5,1]= 6*EI/el**2; ek[5,2]=-12*EI/el**3; ek[5,3]=   0.0; ek[5,4]= 6*EI/el**2; ek[5,5]= 12*EI/el**3
    s=yy/el
    c=xx/el
    tt[0,0]=  c; tt[0,1]=  s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
    tt[1,0]= -s; tt[1,1]=  c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
    tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
    tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]=  c; tt[3,4]=  s; tt[3,5]=0.0
    tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]=  c; tt[4,5]=0.0
    tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
    return ek,tt


def bfvec_grd(gamma,A,x1,y1,x2,y2):
    # Distributed vertical load vector
    bfv=np.zeros(6,dtype=np.float64)
    el=np.sqrt((x2-x1)**2+(y2-y1)**2)
    gkv=-1.0 # downward direction
    bfv[0]=0.0
    bfv[1]=0.0
    bfv[2]=0.5*gamma*A*el*gkv
    bfv[3]=0.0
    bfv[4]=0.0
    bfv[5]=0.5*gamma*A*el*gkv
    return bfv


def main():
    # Main routine
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=3            # Degree of freedom per node
    # data input and input data print
    npoin,nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp=inpdata_grd(fnameR,nod,nfree)
    prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node)
    # assembly of stiffness matrix & load vectors
    ir    =np.zeros(nod*nfree,dtype=np.int)             # Work vector for matrix assembly
    gk    =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64) # Global stiffness matrix
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        aa=ae[2,m]
        ai=ae[3,m]
        aj=ae[4,m]
        gamma=ae[5,m]
        A=aa
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)      # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)     # Stiffness matrix in global coordinate
        bfe  =bfvec_grd(gamma,A,x1,y1,x2,y2) # Body force vector downward direction
        ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
        ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+ck[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of section force
    fsec=np.zeros([nod*nfree,nele],dtype=np.float64) # Section force vector
    dis =np.zeros(nod*nfree,dtype=np.float64) # work vector for section force calculation
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        ai=ae[3,m]
        aj=ae[4,m]
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)
        dis[0]=disg[3*i]; dis[1]=disg[3*i+1]; dis[2]=disg[3*i+2]
        dis[3]=disg[3*j]; dis[4]=disg[3*j+1]; dis[5]=disg[3*j+2]
        fsec[:,ne]=np.dot(ek,np.dot(tt,dis))
    # print out of result
    prout_grd(fnameW,npoin,nele,disg,fsec)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く