Option Explicit On Option Strict On Public Class Form1 Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click Dim sr As System.IO.StreamReader Dim dat As String Dim sbuf() As String Dim delim() As Char = {" "c} Dim fnameL As String = "" 'ファイルリスト Dim fnameR As String = "" '入力ファイル名 Dim fnameW As String = "" '出力ファイル名 '----------------------------------------------------- '入出力データファイルリスト指定 '----------------------------------------------------- OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameL = OpenFileDialog1.FileName '----------------------------------------------------- '読み込み・実行 '----------------------------------------------------- sr = New System.IO.StreamReader(fnameL, System.Text.Encoding.Default) Do Until sr.EndOfStream dat = sr.ReadLine() sbuf = dat.Split(delim) fnameR = System.IO.Path.GetDirectoryName(fnameL) & "\" & sbuf(0) fnameW = System.IO.Path.GetDirectoryName(fnameL) & "\" & sbuf(1) Call main_part(fnameR, fnameW) Loop sr.Close() MessageBox.Show("計算完了", "通知") End Sub Private Sub main_part(ByVal fnameR As String, ByVal fnameW As String) '平面骨組構造解析 Dim i As Integer : Dim j As Integer Dim ne As Integer : Dim n As Integer Dim ia As Integer : Dim ja As Integer Dim sr As System.IO.StreamReader Dim sw As System.IO.StreamWriter Dim dat As String Dim sbuf() As String Dim delim() As Char = {","c} Dim strcom As String Dim NODT As Integer '節点総数 Dim NELT As Integer '要素総数 Dim MATEL As Integer '材料種別数 Dim KOX As Integer 'x方向変位既知節点数 Dim KOY As Integer 'y方向変位既知節点数 Dim KOZ As Integer '回転変位既知節点数 Dim NF As Integer '載荷節点数 Dim nt As Integer '全自由度 Dim mm As Integer '変位既知節点処理語の自由度 Dim kakom(,) As Integer '要素構成接点番号 Dim Em() As Double '要素弾性係数 Dim AA() As Double '要素断面積 Dim AI() As Double '要素断面二次モーメント Dim matno() As Integer '材料種別No Dim x() As Double '節点x座標 Dim y() As Double '節点y座標 Dim nokx() As Integer 'x方向変位既知節点番号 Dim noky() As Integer 'y方向変位既知節点番号 Dim nokz() As Integer '回転変位既知節点番号 Dim df() As Double '全体節点外力増分ベクトル(入力値) Dim ftvec() As Double '全体節点力ベクトル Dim index() As Integer '変位既知節点処理のためのインデックス Dim sme(6, 6) As Double '要素剛性マトリックス Dim smg(6, 6) As Double '要素剛性マトリックス Dim tsm(,) As Double '全体剛性マトリックス Dim tsg(,) As Double '全体剛性マトリックス Dim vec(,) As Double '固有ベクトル Dim dis() As Double '全体節点変位 Dim ev() As Double '固有値 Dim reac() As Double '要素節点反力 Dim hen(6, 6) As Double '座標変換行列 Dim fevec(6) As Double '要素節点力ベクトル Dim wm(,) As Double '材料種別作業用 Dim sttime As Date '計算開始時刻 Dim entime As Date '計算終了時刻 Dim tspan As TimeSpan '計算時間 Dim ib As Integer 'バンド幅 Dim fact0 As Double = 1.0E-30 '0判定 Dim reforce(,) As Double '断面力保存 Dim jdg As Integer '----------------------------------------------------- '時間計測開始 '----------------------------------------------------- sttime = DateTime.Now '***************************************************** 'データ入力 '***************************************************** sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default) '-------------------------------------------------------------------------------------- '書出用コメント入力 '-------------------------------------------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0) '-------------------------------------------------------------------------------------- '節点数,要素数,材料種別数,X方向変位既知節点数,Y方向変位既知節点数,回転変位既知節点数,載荷点数 '-------------------------------------------------------------------------------------- dat = sr.ReadLine : sbuf = dat.Split(delim) NODT = CInt(sbuf(0)) NELT = CInt(sbuf(1)) MATEL = CInt(sbuf(2)) KOX = CInt(sbuf(3)) KOY = CInt(sbuf(4)) KOZ = CInt(sbuf(5)) NF = CInt(sbuf(6)) nt = NODT * 3 '-------------------------------------------------------------------------------------- '配列寸法宣言 '-------------------------------------------------------------------------------------- ReDim kakom(NELT, 2) ReDim Em(NELT) ReDim AA(NELT) ReDim AI(NELT) ReDim matno(NELT) ReDim wm(MATEL, 3) ReDim x(NODT) ReDim y(NODT) ReDim nokx(KOX) ReDim noky(KOY) ReDim nokz(KOZ) ReDim index(nt) ReDim tsm(nt, nt) ReDim tsg(nt, nt) ReDim vec(nt, nt) ReDim df(nt) ReDim ftvec(nt) ReDim dis(nt) ReDim ev(nt) ReDim reac(nt) ReDim reforce(NELT, 6) '------------------------------------------------------------------------------------------ '材料特性(Em,AA,AI) '------------------------------------------------------------------------------------------ For i = 1 To MATEL dat = sr.ReadLine : sbuf = dat.Split(delim) For j = 1 To 3 wm(i, j) = CDbl(sbuf(j - 1)) Next j Next i '------------------------------------------------------------------------------------------ '要素構成節点と材料種別(node-1,node-2,matno) '------------------------------------------------------------------------------------------ For ne = 1 To NELT dat = sr.ReadLine : sbuf = dat.Split(delim) kakom(ne, 1) = CInt(sbuf(0)) kakom(ne, 2) = CInt(sbuf(1)) matno(ne) = CInt(sbuf(2)) For i = 1 To MATEL If i = matno(ne) Then Em(ne) = wm(i, 1) AA(ne) = wm(i, 2) AI(ne) = wm(i, 3) End If Next i Next ne '-------------------------------------------------------------------------------------- '節点座標(x,y) '-------------------------------------------------------------------------------------- For i = 1 To NODT dat = sr.ReadLine : sbuf = dat.Split(delim) x(i) = CDbl(sbuf(0)) y(i) = CDbl(sbuf(1)) Next i '-------------------------------------------------------------------------------------- '変位既知節点番号,既知節点変位量 '-------------------------------------------------------------------------------------- If 0 < KOX Then For i = 1 To KOX dat = sr.ReadLine : sbuf = dat.Split(delim) nokx(i) = CInt(sbuf(0)) Next i End If If 0 < KOY Then For i = 1 To KOY dat = sr.ReadLine : sbuf = dat.Split(delim) noky(i) = CInt(sbuf(0)) Next i End If If 0 < KOZ Then For i = 1 To KOZ dat = sr.ReadLine : sbuf = dat.Split(delim) nokz(i) = CInt(sbuf(0)) Next i End If '-------------------------------------------------------------------------------------- '節点番号,x方向荷重,y方向荷重,モーメント荷重 '-------------------------------------------------------------------------------------- For i = 1 To nt : df(i) = 0.0 : Next i If 0 < NF Then For i = 1 To NF dat = sr.ReadLine : sbuf = dat.Split(delim) n = CInt(sbuf(0)) df(3 * n - 2) = CDbl(sbuf(1)) df(3 * n - 1) = CDbl(sbuf(2)) df(3 * n) = CDbl(sbuf(3)) Next i End If sr.Close() '***************************************************** '入力確認出力 '***************************************************** Call INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOX, KOY, KOZ, NF, _ x, y, df, nokx, noky, nokz, kakom, Em, AA, AI, matno) '***************************************************** 'メイン処理計算 '***************************************************** '----------------------------------------------------- '変位既知節点自由度記憶 '----------------------------------------------------- For i = 1 To nt : index(i) = i : Next i If 0 < KOX Then For i = 1 To KOX : n = 3 * nokx(i) - 2 : index(n) = 0 : Next i End If If 0 < KOY Then For i = 1 To KOY : n = 3 * noky(i) - 1 : index(n) = 0 : Next i End If If 0 < KOZ Then For i = 1 To KOZ : n = 3 * nokz(i) : index(n) = 0 : Next i End If mm = 0 For i = 1 To nt If index(i) > 0 Then mm = mm + 1 : index(mm) = i End If Next i '***************************************************** '初期応力導入 '***************************************************** '初期処理 {df}=(K){dis0} For ne = 1 To NELT For i = 1 To 6 reforce(ne, i) = 0.0 Next i Next ne For i = 1 To nt : For j = 1 To nt : tsm(i, j) = 0.0 : Next j, i For ne = 1 To NELT Call FRAMESM(ne, kakom, x, y, Em, AA, AI, reforce, sme, smg) Call ASMAT(ne, 2, 3, kakom, sme, tsm) Next ne For i = 1 To nt : reac(i) = 0.0 : dis(i) = 0.0 : ftvec(i) = 0.0 : Next i For i = 1 To mm ia = index(i) : reac(i) = df(ia) For j = 1 To mm ja = index(j) : tsm(i, j) = tsm(ia, ja) Next j, i ib = IBAND(mm, tsm, fact0) Call BAND(mm, ib, tsm) jdg = CHBAND10(mm, ib, tsm) '三角行列 If jdg = 1 Then GoTo epart Call CHBAND11(mm, ib, tsm, reac) '前進消去・後退代入 For i = 1 To mm : ia = index(i) : dis(ia) = reac(i) : Next i '内力算定 For ne = 1 To NELT Call CALST(ne, kakom, x, y, dis, Em, AA, AI, fevec, reforce) Call ASVEC(ne, 2, 3, kakom, fevec, ftvec) Next ne Call DATAOUT_ISA(fnameW, NODT, NELT, x, y, dis, df, ftvec, reforce) '***************************************************** '固有値解析 '***************************************************** For i = 1 To nt : For j = 1 To nt : tsm(i, j) = 0.0 : tsg(i, j) = 0.0 : vec(i, j) = 0.0 : Next j, i For ne = 1 To NELT Call FRAMESM(ne, kakom, x, y, Em, AA, AI, reforce, sme, smg) Call ASMAT(ne, 2, 3, kakom, sme, tsm) Call ASMAT(ne, 2, 3, kakom, smg, tsg) Next ne For i = 1 To nt : ev(i) = 0.0 : Next i For i = 1 To mm ia = index(i) For j = 1 To mm ja = index(j) : tsm(i, j) = tsm(ia, ja) : tsg(i, j) = tsg(ia, ja) Next j, i Call GJACOBI(mm, tsm, tsg, ev, vec) For i = 1 To nt dis(i) = 0.0 For j = 1 To nt tsm(i, j) = 0.0 Next j, i For i = 1 To mm ia = index(i) : dis(ia) = ev(i) For j = 1 To mm ja = index(j) : tsm(ja, ia) = vec(j, i) Next j, i DATAOUT_EVA(fnameW, NODT, x, y, dis, tsm) epart: '----------------------------------------------------- '時間計測完了 '----------------------------------------------------- entime = DateTime.Now tspan = entime - sttime sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "NODT=" & NODT.ToString("0") dat = dat & " nt=" & nt.ToString("0") dat = dat & " mm=" & mm.ToString("0") dat = dat & " ib=" & ib.ToString("0") sw.WriteLine(dat) dat = "Calculation Time=" & tspan.TotalSeconds.ToString("0.00sec") sw.WriteLine(dat) dat = entime.ToString() sw.WriteLine(dat) sw.Close() End Sub Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, _ ByVal NODT As Integer, ByVal NELT As Integer, ByVal MATEL As Integer, _ ByVal KOX As Integer, ByVal KOY As Integer, ByVal KOZ As Integer, _ ByVal NF As Integer, ByVal x() As Double, ByRef y() As Double, ByRef df() As Double, _ ByRef nokx() As Integer, ByRef noky() As Integer, ByRef nokz() As Integer, _ ByRef kakom(,) As Integer, ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _ ByRef matno() As Integer) '入力確認出力 Dim sw As System.IO.StreamWriter Dim dat As String Dim i As Integer Dim k As Integer Dim l As Integer Dim m As Integer Dim ne As Integer sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default) dat = strcom : sw.WriteLine(dat) dat = "NODT,NELT,MATEL,KOX,KOY,KOZ,NF" : sw.WriteLine(dat) dat = NODT.ToString("0") dat = dat & "," & NELT.ToString("0") dat = dat & "," & MATEL.ToString("0") dat = dat & "," & KOX.ToString("0") dat = dat & "," & KOY.ToString("0") dat = dat & "," & KOZ.ToString("0") dat = dat & "," & NF.ToString("0") sw.WriteLine(dat) dat = "*input-data" : sw.WriteLine(dat) dat = "node,x,y,fx,fy,fz,x-fix,y-fix,z-fix" : sw.WriteLine(dat) For i = 1 To NODT k = 0 : l = 0 : m = 0 If 0 < KOX Then For j = 1 To KOX If i = nokx(j) Then k = 1 Next j End If If 0 < KOY Then For j = 1 To KOY If i = noky(j) Then l = 1 Next j End If If 0 < KOZ Then For j = 1 To KOZ If i = nokz(j) Then m = 1 Next j End If dat = i.ToString("0") dat = dat & "," & x(i).ToString("E") dat = dat & "," & y(i).ToString("E") dat = dat & "," & df(3 * i - 2).ToString("E") dat = dat & "," & df(3 * i - 1).ToString("E") dat = dat & "," & df(3 * i).ToString("E") dat = dat & "," & k.ToString("0") dat = dat & "," & l.ToString("0") dat = dat & "," & m.ToString("0") sw.WriteLine(dat) Next i dat = "element,node-1,node-2,E,A,I,matno" : sw.WriteLine(dat) For ne = 1 To NELT dat = ne.ToString("0") dat = dat & "," & kakom(ne, 1).ToString("0") dat = dat & "," & kakom(ne, 2).ToString("0") dat = dat & "," & Em(ne).ToString("E") dat = dat & "," & AA(ne).ToString("E") dat = dat & "," & AI(ne).ToString("E") dat = dat & "," & matno(ne).ToString("0") sw.WriteLine(dat) Next ne sw.Close() End Sub Private Sub DATAOUT_ISA(ByVal fnameW As String, ByVal NODT As Integer, ByVal NELT As Integer, _ ByRef x() As Double, ByRef y() As Double, ByRef dis() As Double, _ ByRef df() As Double, ByRef ftvec() As Double, ByRef reforce(,) As Double) '結果出力(Appendモード) Dim i As Integer Dim ne As Integer Dim sw As System.IO.StreamWriter Dim dat As String sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "***** Initial stress analysis *****" : sw.WriteLine(dat) dat = "*Node value" : sw.WriteLine(dat) dat = "node,x-coord,y-coord,x-dir,y-dir,z-dir,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z" : sw.WriteLine(dat) For i = 1 To NODT dat = i.ToString("0") dat = dat & "," & x(i).ToString("E") dat = dat & "," & y(i).ToString("E") dat = dat & "," & dis(3 * i - 2).ToString("E") dat = dat & "," & dis(3 * i - 1).ToString("E") dat = dat & "," & dis(3 * i).ToString("E") dat = dat & "," & df(3 * i - 2).ToString("E") dat = dat & "," & df(3 * i - 1).ToString("E") dat = dat & "," & df(3 * i).ToString("E") dat = dat & "," & ftvec(3 * i - 2).ToString("E") dat = dat & "," & ftvec(3 * i - 1).ToString("E") dat = dat & "," & ftvec(3 * i).ToString("E") sw.WriteLine(dat) Next i dat = "*Stress-resultant" : sw.WriteLine(dat) dat = "element,Ni,Si,Mi,Nj,Sj,Mj" : sw.WriteLine(dat) For ne = 1 To NELT dat = ne.ToString("0") For i = 1 To 6 dat = dat & "," & reforce(ne, i).ToString("E") Next i sw.WriteLine(dat) Next ne sw.Close() End Sub Private Sub DATAOUT_EVA(ByVal fnameW As String, ByVal NODT As Integer, _ ByRef x() As Double, ByRef y() As Double, ByRef EV() As Double, ByRef vec(,) As Double) '結果出力(Appendモード) Dim i As Integer Dim k As Integer Dim nt As Integer Dim sw As System.IO.StreamWriter Dim dat As String nt = NODT * 3 sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "***** Eigen value analysis *****" : sw.WriteLine(dat) i = 0 : k = 0 Do i = i + 1 If 1.0E-30 < Math.Abs(EV(i)) And EV(i) < 0.0 Then k = k + 1 dat = "*Eigen value" : sw.WriteLine(dat) dat = EV(i).ToString("E") : sw.WriteLine(dat) dat = "*Node data" : sw.WriteLine(dat) dat = "i,x,y,vec-x,vec-y,vec-z" : sw.WriteLine(dat) For j = 1 To NODT dat = j.ToString("0") dat = dat & "," & x(j).ToString("E") dat = dat & "," & y(j).ToString("E") dat = dat & "," & vec(3 * j - 2, i).ToString("E") dat = dat & "," & vec(3 * j - 1, i).ToString("E") dat = dat & "," & vec(3 * j, i).ToString("E") sw.WriteLine(dat) Next j End If If 3 <= k Then Exit Do Loop While i < nt sw.Close() End Sub Private Sub FRAMESM(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _ ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _ ByRef reforce(,) As Double, ByRef sme(,) As Double, ByRef smg(,) As Double) Dim i As Integer : Dim j As Integer : Dim k As Integer Dim hen(6, 6) As Double Dim smu(6, 6) As Double Dim smv(6, 6) As Double Dim w(6, 6) As Double Dim AL As Double : Dim P As Double Dim BA As Double : Dim BI As Double Dim xx1 As Double : Dim yy1 As Double : Dim xx2 As Double : Dim yy2 As Double For i = 1 To 6 For j = 1 To 6 smv(i, j) = 0.0 smu(i, j) = 0.0 smg(i, j) = 0.0 sme(i, j) = 0.0 Next j, i '軸力計算 i = kakom(ne, 1) j = kakom(ne, 2) xx1 = x(i) yy1 = y(i) xx2 = x(j) yy2 = y(j) AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)) P = 0.5 * (reforce(ne, 4) - reforce(ne, 1)) '幾何剛性行列 smv(1, 1) = 0.0 smv(2, 1) = 0.0 smv(2, 2) = 6.0 / 5 / AL smv(3, 1) = 0.0 smv(3, 2) = 1.0 / 10.0 smv(3, 3) = 2.0 * AL / 15.0 smv(4, 1) = 0.0 smv(4, 2) = 0.0 smv(4, 3) = 0.0 smv(4, 4) = 0.0 smv(5, 1) = 0.0 smv(5, 2) = -6.0 / 5.0 / AL smv(5, 3) = -1.0 / 10.0 smv(5, 4) = 0.0 smv(5, 5) = 6.0 / 5.0 / AL smv(6, 1) = 0.0 smv(6, 2) = 1.0 / 10.0 smv(6, 3) = -AL / 30.0 smv(6, 4) = 0.0 smv(6, 5) = -1.0 / 10.0 smv(6, 6) = 2.0 * AL / 15.0 For i = 1 To 6 For j = i To 6 smv(i, j) = smv(j, i) Next j, i BA = AA(ne) BI = AI(ne) smu(1, 1) = 1.0 / AL smu(2, 1) = 0.0 smu(2, 2) = 12.0 * BI / BA / AL / AL / AL smu(3, 1) = 0.0 smu(3, 2) = 6 * BI / BA / AL / AL smu(3, 3) = 4 * BI / BA / AL smu(4, 1) = -1.0 / AL smu(4, 2) = 0.0 smu(4, 3) = 0.0 smu(4, 4) = 1.0 / AL smu(5, 1) = 0.0 smu(5, 2) = -12.0 * BI / BA / AL / AL / AL smu(5, 3) = -6.0 * BI / BA / AL / AL smu(5, 4) = 0.0 smu(5, 5) = 12.0 * BI / BA / AL / AL / AL smu(6, 1) = 0.0 smu(6, 2) = 6.0 * BI / BA / AL / AL smu(6, 3) = 2.0 * BI / BA / AL smu(6, 4) = 0.0 smu(6, 5) = 6.0 * BI / BA / AL / AL smu(6, 6) = 4.0 * BI / BA / AL For i = 1 To 6 For j = i To 6 smu(i, j) = smu(j, i) Next j, i For i = 1 To 6 For j = 1 To 6 smg(i, j) = smg(i, j) + P * smv(i, j) + P * smu(i, j) Next j, i Call HENMX(ne, xx1, yy1, xx2, yy2, hen) '(w)=(T)'*(k) For i = 1 To 6 For j = 1 To 6 w(i, j) = 0.0 For k = 1 To 6 w(i, j) = w(i, j) + hen(k, i) * smg(k, j) Next k, j, i '(w)*(T)=(T)'(k,T):(sm)は全体座標系での要素剛性行列 For i = 1 To 6 For j = 1 To 6 smg(i, j) = 0.0 For k = 1 To 6 smg(i, j) = smg(i, j) + w(i, k) * hen(k, j) Next k, j, i '線形弾性剛性行列 Call ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme) '(w)=(T)'*(k) For i = 1 To 6 For j = 1 To 6 w(i, j) = 0.0 For k = 1 To 6 w(i, j) = w(i, j) + hen(k, i) * sme(k, j) Next k, j, i '(w)*(T)=(T)'(k,T):(sm)は全体座標系での要素剛性行列 For i = 1 To 6 For j = 1 To 6 sme(i, j) = 0.0 For k = 1 To 6 sme(i, j) = sme(i, j) + w(i, k) * hen(k, j) Next k, j, i End Sub Private Sub ELMSM(ByVal ne As Integer, _ ByVal xx1 As Double, ByVal yy1 As Double, ByVal xx2 As Double, ByVal yy2 As Double, _ ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, ByRef sme(,) As Double) '要素剛性行列 Dim i As Integer : Dim j As Integer Dim AL As Double : Dim EA As Double : Dim EI As Double AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)) EA = Em(ne) * AA(ne) EI = Em(ne) * AI(ne) For i = 1 To 6 For j = 1 To 6 sme(i, j) = 0.0 Next j, i sme(1, 1) = EA / AL sme(2, 2) = 12.0 * EI / (AL * AL * AL) sme(3, 2) = 6.0 * EI / (AL * AL) sme(3, 3) = 4.0 * EI / AL sme(4, 1) = -EA / AL sme(4, 4) = sme(1, 1) sme(5, 2) = -sme(2, 2) sme(5, 3) = -sme(3, 2) sme(5, 5) = sme(2, 2) sme(6, 2) = sme(3, 2) sme(6, 3) = 2.0 * EI / AL sme(6, 5) = -sme(3, 2) sme(6, 6) = sme(3, 3) '対称項セット For i = 1 To 6 For j = 1 To i sme(j, i) = sme(i, j) Next j, i End Sub Private Sub CALST(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _ ByRef dis() As Double, ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _ ByRef fevec() As Double, ByRef reforce(,) As Double) '断面力算定 Dim i As Integer : Dim k As Integer Dim ii As Integer : Dim jj As Integer : Dim ia As Integer Dim hen(6, 6) As Double Dim sme(6, 6) As Double Dim wd(6) As Double : Dim wde(6) As Double : Dim ws(6) As Double Dim xx1 As Double : Dim yy1 As Double : Dim xx2 As Double : Dim yy2 As Double ii = kakom(ne, 1) jj = kakom(ne, 2) xx1 = x(ii) yy1 = y(ii) xx2 = x(jj) yy2 = y(jj) Call HENMX(ne, xx1, yy1, xx2, yy2, hen) Call ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme) For i = 1 To 2 ia = kakom(ne, i) wd(3 * i - 2) = dis(3 * ia - 2) wd(3 * i - 1) = dis(3 * ia - 1) wd(3 * i) = dis(3 * ia) Next i For i = 1 To 6 wde(i) = 0.0 For k = 1 To 6 wde(i) = wde(i) + hen(i, k) * wd(k) Next k, i For i = 1 To 6 ws(i) = 0.0 For k = 1 To 6 ws(i) = ws(i) + sme(i, k) * wde(k) Next k, i '断面力保存 For i = 1 To 6 reforce(ne, i) = reforce(ne, i) + ws(i) Next i For i = 1 To 6 fevec(i) = 0.0 For k = 1 To 6 fevec(i) = fevec(i) + hen(k, i) * reforce(ne, k) Next k, i End Sub Private Sub HENMX(ByVal ne As Integer, ByRef xx1 As Double, ByVal yy1 As Double, ByVal xx2 As Double, ByVal yy2 As Double, ByRef hen(,) As Double) '座標変換行列 Dim i As Integer Dim j As Integer Dim AL As Double Dim cs As Double Dim sn As Double AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)) cs = (xx2 - xx1) / AL sn = (yy2 - yy1) / AL For i = 1 To 6 For j = 1 To 6 hen(i, j) = 0.0 Next j, i hen(1, 1) = cs hen(1, 2) = sn hen(2, 1) = -sn hen(2, 2) = cs hen(3, 3) = 1.0 For i = 1 To 3 For j = 1 To 3 hen(i + 3, j + 3) = hen(i, j) Next j, i End Sub Private Sub ASMAT(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _ ByRef kakom(,) As Integer, ByRef sm(,) As Double, ByRef tsm(,) As Double) '剛性行列組立(nod:1要素節点数,nhen:1節点の自由度) Dim i As Integer : Dim j As Integer : Dim k As Integer : Dim l As Integer Dim ki As Integer : Dim kj As Integer : Dim ks As Integer : Dim js As Integer Dim kik As Integer : Dim kjl As Integer : Dim ksk As Integer : Dim jsl As Integer For i = 1 To nod For j = 1 To nod ki = (kakom(ne, i) - 1) * nhen kj = (kakom(ne, j) - 1) * nhen ks = (i - 1) * nhen js = (j - 1) * nhen For k = 1 To nhen For l = 1 To nhen kik = ki + k : kjl = kj + l : ksk = ks + k : jsl = js + l tsm(kik, kjl) = tsm(kik, kjl) + sm(ksk, jsl) Next l, k Next j, i End Sub Private Sub ASVEC(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _ ByRef kakom(,) As Integer, ByRef fevec() As Double, ByRef ftvec() As Double) '荷重ベクトル組立(nod:1要素節点数,nhen:1節点の自由度) Dim i As Integer : Dim k As Integer Dim ki As Integer : Dim ks As Integer : Dim kik As Integer : Dim ksk As Integer For i = 1 To nod ki = (kakom(ne, i) - 1) * nhen ks = (i - 1) * nhen For k = 1 To nhen kik = ki + k : ksk = ks + k ftvec(kik) = ftvec(kik) + fevec(ksk) Next k Next i End Sub Private Function CHBAND10(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double) As Integer Dim i As Integer Dim j As Integer Dim k As Integer Dim mm As Integer Dim mn As Integer Dim z As Double array(1, 1) = Math.Sqrt(array(1, 1)) For j = 2 To m1 array(1, j) = array(1, j) / array(1, 1) Next j For i = 2 To n1 If n1 - i + 1 > m1 Then mm = m1 Else mm = n1 - i + 1 End If For j = 1 To mm z = array(i, j) If j <> m1 Then If m1 - j + 1 < i Then mn = m1 - j + 1 Else mn = i End If For k = 2 To mn z = z - array(i - k + 1, k) * array(i - k + 1, j + k - 1) Next k If j > 1 Then array(i, j) = z / array(i, 1) Else If z < 0.0 Then Return 1 array(i, 1) = Math.Sqrt(z) End If Else array(i, j) = z / array(i, 1) End If Next j, i Return 0 End Function Private Sub CHBAND11(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double, ByRef vector() As Double) Dim i As Integer Dim j As Integer Dim mm As Integer Dim z As Double vector(1) = vector(1) / array(1, 1) For i = 2 To n1 z = vector(i) If i > m1 Then mm = m1 Else mm = i End If For j = 2 To mm z = z - array(i - j + 1, j) * vector(i - j + 1) Next j vector(i) = z / array(i, 1) Next i vector(n1) = vector(n1) / array(n1, 1) For i = n1 - 1 To 1 Step -1 z = vector(i) For j = 2 To m1 If i + j - 1 > n1 Then Exit For z = z - array(i, j) * vector(i + j - 1) Next j vector(i) = z / array(i, 1) Next i End Sub Private Function IBAND(ByVal mm As Integer, ByRef array(,) As Double, ByVal fact0 As Double) As Integer Dim i As Integer Dim j As Integer '----------------------------------------------------- 'バンド幅計算 '----------------------------------------------------- IBAND = 1 For i = 1 To mm - 1 For j = mm To i + 1 Step -1 If fact0 <= Math.Abs(array(i, j)) Then Exit For Next j If IBAND <= j - i + 1 Then IBAND = j - i + 1 Next i If mm < IBAND Then IBAND = mm End Function Private Sub BAND(ByVal mm As Integer, ByVal ib As Integer, ByRef array(,) As Double) Dim i As Integer Dim j As Integer Dim mn As Integer Dim work() As Double ReDim work(mm) '----------------------------------------------------- 'フルマトリックスのバンド化記憶 '----------------------------------------------------- For i = 2 To mm If i + ib - 1 < mm Then mn = ib Else mn = mm - i + 1 End If For j = 1 To ib work(j) = 0.0 Next j For k = 1 To mn work(k) = array(i, i + k - 1) Next k For j = 1 To ib array(i, j) = work(j) Next j Next i End Sub Private Sub IERROR1(ByVal fnameW As String) '結果出力(Appendモード) Dim sw As System.IO.StreamWriter Dim dat As String sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "area of element is less or equal to 0" : sw.WriteLine(dat) sw.Close() MessageBox.Show("要素面積<=0", "異常終了通知") Me.Close() End Sub Private Sub IERROR2(ByVal fnameW As String) '結果出力(Appendモード) Dim sw As System.IO.StreamWriter Dim dat As String sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default) dat = "diagonal value is less or equal to 0" : sw.WriteLine(dat) sw.Close() MessageBox.Show("対角項<=0", "異常終了通知") Me.Close() End Sub Private Sub GJACOBI(ByVal nnf As Integer, ByRef a(,) As Double, ByRef b(,) As Double, ByRef ev() As Double, ByRef vec(,) As Double) '************************************ '一般化JACOBI法による固有値計算 '(パソコン構造力学,宮澤健二著) '************************************ Dim i As Integer Dim j As Integer Dim k As Integer Dim fmk As Double Dim fmg As Double Dim fact As Double = 1.0E-20 Dim epsk As Double Dim epsg As Double Dim kmax As Integer Dim ip As Integer Dim iq As Integer Dim amax As Double Dim dam As Double Dim de1 As Double Dim de2 As Double Dim de3 As Double Dim dd As Double Dim ss As Double Dim x As Double Dim gg As Double Dim ww As Double Dim im As Integer fmk = 0.0 fmg = 0.0 For i = 1 To nnf For j = 1 To nnf If Math.Abs(a(i, j)) > fmk Then fmk = Math.Abs(a(i, j)) If Math.Abs(b(i, j)) > fmg Then fmg = Math.Abs(b(i, j)) Next j Next i epsk = fmk * fact epsg = fmg * fact For i = 1 To nnf For j = 1 To nnf vec(i, j) = 0.0 Next j vec(i, i) = 1.0 Next i kmax = 5 * nnf * nnf For k = 1 To kmax ip = 1 iq = 2 amax = a(1, 2) ^ 2 / (a(1, 1) ^ 2 + a(2, 2) ^ 2) For i = 1 To nnf - 1 For j = i + 1 To nnf dam = a(i, j) ^ 2 / (a(i, i) ^ 2 + a(j, j) ^ 2) If amax < dam Then ip = i : iq = j : amax = dam End If Next j Next i For i = 1 To nnf - 1 For j = i + 1 To nnf dam = b(i, j) ^ 2 / (b(i, i) ^ 2 + b(j, j) ^ 2) If amax < dam Then ip = i : iq = j : amax = dam End If Next j Next i de1 = a(ip, ip) * b(ip, iq) - a(ip, iq) * b(ip, ip) de2 = a(iq, iq) * b(ip, iq) - a(ip, iq) * b(iq, iq) de3 = a(ip, ip) * b(iq, iq) - a(iq, iq) * b(ip, ip) dd = de3 * de3 + 4.0 * de1 * de2 ss = Math.Sqrt(dd) If de3 < 0 Then x = 0.5 * (de3 - ss) Else x = 0.5 * (de3 + ss) End If If x = 0.0 Then dam = 0.0 : gg = -a(ip, iq) / a(iq, iq) Else dam = de2 / x : gg = -de1 / x End If For j = 1 To nnf ww = a(ip, j) a(ip, j) = a(ip, j) + gg * a(iq, j) a(iq, j) = a(iq, j) + dam * ww ww = b(ip, j) b(ip, j) = b(ip, j) + gg * b(iq, j) b(iq, j) = b(iq, j) + dam * ww Next j For i = 1 To nnf ww = a(i, ip) a(i, ip) = a(i, ip) + gg * a(i, iq) a(i, iq) = a(i, iq) + dam * ww ww = b(i, ip) b(i, ip) = b(i, ip) + gg * b(i, iq) b(i, iq) = b(i, iq) + dam * ww ww = vec(i, ip) vec(i, ip) = vec(i, ip) + gg * vec(i, iq) vec(i, iq) = vec(i, iq) + dam * ww Next i If amax < epsk And amax < epsg Then Exit For Next k For i = 1 To nnf ev(i) = a(i, i) / b(i, i) Next i '並び替え 'ev(k):最大固有値k=1,最小固有値k=nnf 'vec(i,k):k番目の固有値に対応する固有ベクトル For k = 1 To nnf - 1 amax = -1.0E+30 For i = k To nnf If ev(i) >= amax Then im = i : amax = ev(i) End If Next i ww = ev(k) : ev(k) = ev(im) : ev(im) = ww For i = 1 To nnf ww = vec(i, k) vec(i, k) = vec(i, im) vec(i, im) = ww Next i Next k End Sub End Class