Option Explicit On Option Strict On Public Class Form2 Private KBD As Integer = 500 '最大境界数 Private KTJ As Integer = 10000 '最大節点数 Private KCM As Integer = 100 '節点に集まる最大要素数 Private KTE As Integer Private LTE As Integer = 100 '多角形の最大辺数 Private NODT As Integer '節点数(入力値) Private pi As Double = Math.PI Private Sub Form2_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load Call MAIN() Me.Close() End Sub Private Sub MAIN() Dim ibex() As Integer '個々の外部境界上の節点数 Dim ibin() As Integer '個々の内部境界上の節点数 Dim ibno(,) As Integer '外部境界上の節点番号 Dim index() As Integer '作業領域 Dim px() As Double '節点のx座標 Dim py() As Double '節点のy座標 Dim mtj(,) As Integer '要素-節点関係 Dim jac(,) As Integer '要素-隣接要素関係 Dim idm() As Integer '要素-領域関係 ReDim ibex(KBD) ReDim ibin(KBD) ReDim ibno(KBD, KTJ) ReDim index(KBD + 1) ReDim px(KTJ + 3) ReDim py(KTJ + 3) ReDim mtj(2 * KTJ + 1, 3) ReDim jac(2 * KTJ + 1, 3) ReDim idm(2 * KTJ + 1) Dim NEX As Integer '外部境界数 Dim NIN As Integer '内部境界数 Dim NOB As Integer '境界上の節点総数 Dim NIB As Integer '境界内部の節点総数 Dim NODE As Integer '節点総数 Dim NELM As Integer '要素総数 Call INPUT(NEX, NIN, ibex, ibin, ibno, NOB, NIB, px, py) Call MODEL(NEX, NIN, ibex, ibin, ibno, NOB, NIB, index, NODE, px, py, NELM, mtj, jac, idm) Call DATA(NEX, NODE, NELM, mtj, jac, idm, px, py) End Sub Private Sub INPUT(ByRef NEX As Integer, ByRef NIN As Integer, ByRef ibex() As Integer, ByRef ibin() As Integer, _ ByRef ibno(,) As Integer, ByRef NOB As Integer, ByRef NIB As Integer, _ ByRef px() As Double, ByRef py() As Double) Dim sr As System.IO.StreamReader Dim dat As String Dim sbuf() As String Dim delim() As Char = {","c} Dim fname1 As String Dim i As Integer Dim j As Integer fname1 = My.Forms.Form1.Label1.Text sr = New System.IO.StreamReader(fname1, System.Text.Encoding.Default) dat = sr.ReadLine() : sbuf = dat.Split(delim) NODT = CInt(sbuf(0)) dat = sr.ReadLine() : sbuf = dat.Split(delim) NEX = CInt(sbuf(0)) NIN = CInt(sbuf(1)) If 0 < NEX Then For i = 1 To NEX dat = sr.ReadLine() : sbuf = dat.Split(delim) ibex(i) = CInt(sbuf(0)) Next i End If If 0 < NIN Then For i = 1 To NIN dat = sr.ReadLine() : sbuf = dat.Split(delim) ibin(i) = CInt(sbuf(0)) Next i End If For i = 1 To NEX dat = sr.ReadLine() : sbuf = dat.Split(delim) For j = 1 To ibex(i) ibno(i, j) = CInt(sbuf(j - 1)) Next j, i dat = sr.ReadLine() : sbuf = dat.Split(delim) NOB = CInt(sbuf(0)) NIB = CInt(sbuf(1)) For i = 1 To NOB + NIB dat = sr.ReadLine() : sbuf = dat.Split(delim) px(i) = CDbl(sbuf(0)) py(i) = CDbl(sbuf(1)) Next i sr.Close() End Sub Private Sub MODEL(ByRef NEX As Integer, ByRef NIN As Integer, ByRef ibex() As Integer, ByRef ibin() As Integer, _ ByRef ibno(,) As Integer, ByRef NOB As Integer, ByRef NIB As Integer, _ ByRef index() As Integer, ByRef NODE As Integer, ByRef px() As Double, ByRef py() As Double, _ ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef idm() As Integer) Dim i As Integer Dim j As Integer Dim ntp As Integer Dim xmin As Double Dim xmax As Double Dim ymin As Double Dim ymax As Double Dim rax As Double Dim ray As Double Dim dmax As Double Dim ia As Integer Dim ib As Integer Dim ic As Integer Dim jnb() As Integer Dim nei(,) As Integer Dim ifix() As Integer Dim angl() As Double Dim nedg() As Integer Dim kstack() As Integer Dim list() As Integer Dim istack() As Integer Dim kv() As Integer Dim ibr() As Integer Dim iadres() As Integer Dim map() As Integer ReDim jnb(KTJ) ReDim nei(KTJ, KCM) ReDim ifix(KTJ) ReDim angl(2 * KTJ + 1) ReDim nedg(2 * KTJ + 1) ReDim kstack(2 * KTJ + 1) ReDim list(KTJ) ReDim istack(KTJ) ReDim kv(2 * KTJ + 1) ReDim ibr(2 * KTJ + 1) ReDim iadres(KTJ + 3) ReDim map(2 * KTJ + 1) '初期化 For i = 1 To KBD + 1 index(i) = 0 Next i KTE = 2 * KTJ + 1 NODE = 0 NELM = 0 For i = 1 To KTE idm(i) = 0 For j = 1 To 3 mtj(i, j) = 0 jac(i, j) = 0 Next j, i For i = 1 To KTJ jnb(i) = 0 For j = 1 To KCM nei(i, j) = 0 Next j, i For i = 1 To KTJ ifix(i) = 0 Next i ntp = NOB + NIB xmin = px(1) xmax = xmin ymin = py(1) ymax = ymin For i = 2 To ntp If px(i) < xmin Then xmin = px(i) If xmax < px(i) Then xmax = px(i) If py(i) < ymin Then ymin = py(i) If ymax < py(i) Then ymax = py(i) Next i rax = xmax - xmin ray = ymax - ymin dmax = rax If dmax < ray Then dmax = ray For i = 1 To ntp px(i) = (px(i) - xmin) / dmax py(i) = (py(i) - ymin) / dmax Next i NELM = 1 ia = KTJ + 1 ib = KTJ + 2 ic = KTJ + 3 mtj(1, 1) = ia mtj(1, 2) = ib mtj(1, 3) = ic jac(1, 1) = 0 jac(1, 2) = 0 jac(1, 3) = 0 px(ia) = -1.23 py(ia) = -0.5 px(ib) = 2.23 py(ib) = -0.5 px(ic) = 0.5 py(ic) = 2.5 Call ROUGH(NEX, NIN, ibex, ibin, ibno, ntp, NIB, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, _ list, iadres, istack, kv, ibr, map) For i = 1 To NODE ifix(i) = 1 Next i Call FINE(NODE, px, py, jnb, nei, ifix, NELM, mtj, jac, idm, istack, angl, nedg, kstack, map) If NELM <> 2 * NODE + 1 Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB MODEL", "Error終了") Me.Close() End If Call LAPLAS(NODE, ifix, mtj, px, py, jnb, nei) Call REMOVE(NEX, index, NELM, mtj, jac, idm) Call CHECK(NELM, mtj, jac) For i = 1 To NODE px(i) = px(i) * dmax + xmin py(i) = py(i) * dmax + ymin Next i End Sub Private Sub ROUGH(ByRef NEX As Integer, ByRef NIN As Integer, ByRef ibex() As Integer, ByRef ibin() As Integer, _ ByRef ibno(,) As Integer, ByRef ntp As Integer, ByRef NIB As Integer, _ ByRef NODE As Integer, ByRef px() As Double, ByRef py() As Double, _ ByRef jnb() As Integer, ByRef nei(,) As Integer, _ ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef idm() As Integer, _ ByRef list() As Integer, ByRef iadres() As Integer, ByRef istack() As Integer, ByRef kv() As Integer, _ ByRef ibr() As Integer, ByRef map() As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim nb As Integer Dim np As Integer Dim ma As Integer Dim mb As Integer Dim mc As Integer Dim ms As Integer Dim mp As Integer Dim js As Integer Dim jz As Integer Dim xs As Double Dim ys As Double Dim loc As Integer Dim xa As Double Dim ya As Double Dim it As Integer For i = 1 To NEX jz = 0 nb = i np = ibex(nb) For j = 1 To np list(j) = ibno(nb, j) Next j Call BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE) For k = 1 To NELM If idm(k) = jz Then ma = iadres(mtj(k, 1)) mb = iadres(mtj(k, 2)) mc = iadres(mtj(k, 3)) ms = ma * mb * mc mp = (mb - ma) * (mc - mb) * (ma - mc) If (ms <> 0) And (0 < mp) Then idm(k) = nb End If End If Next k Next i For i = 1 To NIN nb = i np = ibin(nb) For j = 1 To np list(j) = NODE + j Next j js = list(1) xs = px(js) ys = py(js) Call LOCATE(xs, ys, px, py, mtj, jac, NELM, loc) jz = idm(loc) Call BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE) For k = 1 To NELM If idm(k) = jz Then ma = iadres(mtj(k, 1)) mb = iadres(mtj(k, 2)) mc = iadres(mtj(k, 3)) ms = ma * mb * mc mp = (mb - ma) * (mc - mb) * (ma - mc) If (ms <> 0) And (mp < 0) Then idm(k) = 0 End If End If Next k Next i For i = 1 To KTJ + 3 iadres(i) = 0 Next i For i = 1 To NIB NODE = NODE + 1 xa = px(NODE) ya = py(NODE) Call LOCATE(xa, ya, px, py, mtj, jac, NELM, it) jz = idm(it) Call DELAUN(jz, NODE, NODE, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack) Next i If NODE <> ntp Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB ROUGH", "Error終了") Me.Close() End If End Sub Private Sub BOUGEN(ByRef jz As Integer, ByRef np As Integer, ByRef list() As Integer, ByRef ntp As Integer, _ ByRef px() As Double, ByRef py() As Double, ByRef jnb() As Integer, ByRef nei(,) As Integer, _ ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef idm() As Integer, _ ByRef iadres() As Integer, ByRef istack() As Integer, ByRef kv() As Integer, ByRef ibr() As Integer, _ ByRef map() As Integer, ByRef NODE As Integer) Dim i As Integer Dim inp As Integer Dim js As Integer Dim ip As Integer Dim iq As Integer Dim iv As Integer inp = 0 For i = 1 To KTJ + 3 iadres(i) = 0 Next i For i = 1 To np iadres(list(i)) = i Next i js = list(1) If NODE < js Then inp = inp + 1 Call DELAUN(jz, js, js, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack) End If For i = 1 To np ip = list((i Mod np) + 1) iq = list(i) If (NODE < ip) And (i <> np) Then inp = inp + 1 Call DELAUN(jz, ip, ip, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack) End If For j = 1 To jnb(iq) For k = 1 To 3 If mtj(nei(iq, j), k) = ip Then GoTo bougen30 Next k Next j Call SEARCH(jz, ip, iq, jnb, nei, NELM, mtj, jac, idm, istack, iv, kv, iadres, ibr) Call POLY(iq, ip, iv, kv, px, py, NELM, mtj, jac, jnb, nei, map) bougen30: Next i NODE = NODE + inp End Sub Private Sub DELAUN(ByRef jz As Integer, ByRef js As Integer, ByRef jg As Integer, ByRef ntp As Integer, ByRef px() As Double, ByRef py() As Double, _ ByRef jnb() As Integer, ByRef nei(,) As Integer, ByRef NELM As Integer, ByRef mtj(,) As Integer, _ ByRef jac(,) As Integer, ByRef idm() As Integer, ByRef iadres() As Integer, ByRef istack() As Integer) Dim itop As Integer Dim maxstk As Integer Dim ip As Integer Dim it As Integer Dim xp As Double Dim yp As Double Dim ia As Integer Dim ib As Integer Dim ic As Integer Dim iv1 As Integer Dim iv2 As Integer Dim iv3 As Integer Dim idf As Integer Dim ms As Integer Dim iedge As Integer Dim il As Integer Dim ir As Integer Dim iera As Integer Dim ierb As Integer Dim ierl As Integer Dim iswap As Integer itop = 0 maxstk = ntp For i = js To jg ip = i xp = px(ip) yp = py(ip) Call LOCATE(xp, yp, px, py, mtj, jac, NELM, it) If idm(it) <> jz Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB DELAUN", "Error終了") Me.Close() End If ia = jac(it, 1) ib = jac(it, 2) ic = jac(it, 3) iv1 = mtj(it, 1) iv2 = mtj(it, 2) iv3 = mtj(it, 3) mtj(it, 1) = ip mtj(it, 2) = iv1 mtj(it, 3) = iv2 jac(it, 1) = NELM + 2 jac(it, 2) = ia jac(it, 3) = NELM + 1 NELM = NELM + 1 idm(NELM) = jz mtj(NELM, 1) = ip mtj(NELM, 2) = iv2 mtj(NELM, 3) = iv3 jac(NELM, 1) = it jac(NELM, 2) = ib jac(NELM, 3) = NELM + 1 NELM = NELM + 1 idm(NELM) = jz mtj(NELM, 1) = ip mtj(NELM, 2) = iv3 mtj(NELM, 3) = iv1 jac(NELM, 1) = NELM - 1 jac(NELM, 2) = ic jac(NELM, 3) = it Call INCR(iv1, NELM, jnb, nei) Call INCR(iv2, NELM - 1, jnb, nei) If iv3 <= KTJ Then nei(iv3, NEIBOR(iv3, it, jnb, nei)) = NELM - 1 End If Call INCR(iv3, NELM, jnb, nei) jnb(ip) = 3 nei(ip, 1) = it nei(ip, 2) = NELM - 1 nei(ip, 3) = NELM If ia <> 0 Then ms = iadres(iv1) * iadres(iv2) idf = Math.Abs(iadres(iv1) - iadres(iv2)) If (idm(ia) = jz) And ((ms = 0) Or (idf <> 1)) Then itop = itop + 1 istack(itop) = IPUSH(it, maxstk, itop, istack) End If End If If ib <> 0 Then Call EDGE(ib, it, jac, iedge) jac(ib, iedge) = NELM - 1 ms = iadres(iv2) * iadres(iv3) idf = Math.Abs(iadres(iv2) - iadres(iv3)) If (idm(ib) = jz) And ((ms = 0) Or (idf <> 1)) Then itop = itop + 1 istack(itop) = IPUSH(NELM - 1, maxstk, itop, istack) End If End If If ic <> 0 Then Call EDGE(ic, it, jac, iedge) jac(ic, iedge) = NELM ms = iadres(iv3) * iadres(iv1) idf = Math.Abs(iadres(iv3) - iadres(iv1)) If (idm(ic) = jz) And ((ms = 0) Or (idf <> 1)) Then itop = itop + 1 istack(itop) = IPUSH(NELM, maxstk, itop, istack) End If End If delaun50: If 0 < itop Then il = istack(itop) itop = itop - 1 ir = jac(il, 2) Call EDGE(ir, il, jac, ierl) iera = (ierl Mod 3) + 1 ierb = (iera Mod 3) + 1 iv1 = mtj(ir, ierl) iv2 = mtj(ir, iera) iv3 = mtj(ir, ierb) Call SWAP(px(iv1), py(iv1), px(iv2), py(iv2), px(iv3), py(iv3), xp, yp, iswap) If iswap = 1 Then ia = jac(ir, iera) ib = jac(ir, ierb) ic = jac(il, 3) mtj(il, 3) = iv3 jac(il, 2) = ia jac(il, 3) = ir mtj(ir, 1) = ip mtj(ir, 2) = iv3 mtj(ir, 3) = iv1 jac(ir, 1) = il jac(ir, 2) = ib jac(ir, 3) = ic Call DECR(iv1, il, jnb, nei) Call DECR(iv2, ir, jnb, nei) Call INCR(ip, ir, jnb, nei) Call INCR(iv3, il, jnb, nei) If ia <> 0 Then Call EDGE(ia, ir, jac, iedge) jac(ia, iedge) = il ms = iadres(iv2) * iadres(iv3) idf = Math.Abs(iadres(iv2) - iadres(iv3)) If (idm(ia) = jz) And ((ms = 0) Or (idf <> 1)) Then itop = itop + 1 istack(itop) = IPUSH(il, maxstk, itop, istack) End If End If If ib <> 0 Then ms = iadres(iv3) * iadres(iv1) idf = Math.Abs(iadres(iv3) - iadres(iv1)) If (idm(ib) = jz) And ((ms = 0) Or (idf <> 1)) Then itop = itop + 1 istack(itop) = IPUSH(ir, maxstk, itop, istack) End If End If If ic <> 0 Then Call EDGE(ic, il, jac, iedge) jac(ic, iedge) = ir End If End If GoTo delaun50 End If Next i End Sub Private Sub SEARCH(ByRef iz As Integer, ByRef ip As Integer, ByRef iq As Integer, ByRef jnb() As Integer, _ ByRef nei(,) As Integer, ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, _ ByRef idm() As Integer, ByRef istack() As Integer, ByRef iv As Integer, _ ByRef kv() As Integer, ByRef iadres() As Integer, ByRef ibr() As Integer) Dim i As Integer Dim j As Integer Dim ia As Integer Dim ib As Integer Dim ja As Integer Dim jb As Integer Dim mstk As Integer Dim nbr As Integer Dim idf As Integer Dim ms As Integer Dim ielm As Integer Dim jelm As Integer Dim kelm As Integer Dim min As Integer Dim jr As Integer iv = 0 mstk = 0 nbr = 0 For i = 1 To NELM kv(i) = 0 ibr(i) = 0 Next i For i = 1 To KTJ istack(i) = 0 Next i For i = 1 To jnb(iq) nbr = nbr + 1 ibr(nei(iq, i)) = nbr Next i For i = 1 To jnb(iq) ielm = nei(iq, i) j = IVERT(ielm, iq, mtj) ja = (j Mod 3) + 1 jb = (ja Mod 3) + 1 ia = mtj(ielm, ja) ib = mtj(ielm, jb) idf = Math.Abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) If (idf = 1) And (ms <> 0) Then GoTo search40 jelm = jac(ielm, ja) If jelm = 0 Then GoTo search40 If idm(jelm) <> iz Then GoTo search40 nbr = nbr + 1 ibr(jelm) = nbr If mtj(jelm, 1) = ip Then GoTo search80 If mtj(jelm, 2) = ip Then GoTo search80 If mtj(jelm, 3) = ip Then GoTo search80 mstk = mstk + 1 istack(mstk) = jelm search40: Next i search50: If mstk = 0 Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB SEARCH", "Error終了") Me.Close() End If ielm = istack(1) mstk = mstk - 1 For i = 1 To mstk istack(i) = istack(i + 1) Next i istack(mstk + 1) = 0 For j = 1 To 3 ja = (j Mod 3) + 1 ia = mtj(ielm, j) ib = mtj(ielm, ja) idf = Math.Abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) If (idf = 1) And (ms <> 0) Then GoTo search70 jelm = jac(ielm, j) If jelm = 0 Then GoTo search70 If ibr(jelm) <> 0 Then GoTo search70 If idm(jelm) <> iz Then GoTo search70 nbr = nbr + 1 ibr(jelm) = nbr If mtj(jelm, 1) = ip Then GoTo search80 If mtj(jelm, 2) = ip Then GoTo search80 If mtj(jelm, 3) = ip Then GoTo search80 mstk = mstk + 1 istack(mstk) = jelm search70: Next j GoTo search50 search80: iv = iv + 1 kv(iv) = jelm If ibr(jelm) <= jnb(iq) Then GoTo search100 min = KTE + 1 For j = 1 To 3 jr = jac(jelm, j) If jr = 0 Then GoTo search90 If ibr(jr) = 0 Then GoTo search90 ja = (j Mod 3) + 1 ia = mtj(jelm, j) ib = mtj(jelm, ja) idf = Math.Abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) If (idf = 1) And (ms <> 0) Then GoTo search90 If ibr(jr) < min Then kelm = jr min = ibr(jr) End If search90: Next j If min = KTE + 1 Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB SEARCH", "Error終了") Me.Close() End If jelm = kelm GoTo search80 search100: End Sub Private Sub POLY(ByRef iq As Integer, ByRef ip As Integer, ByRef iv As Integer, ByRef kv() As Integer, _ ByRef px() As Double, ByRef py() As Double, ByRef NELM As Integer, ByRef mtj(,) As Integer, _ ByRef jac(,) As Integer, ByRef jnb() As Integer, ByRef nei(,) As Integer, ByRef map() As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim l As Integer Dim ia As Integer Dim ja As Integer Dim ir As Integer Dim jr As Integer Dim iv1 As Integer Dim iv2 As Integer Dim iedge As Integer Dim npa As Integer Dim npb As Integer Dim nta As Integer Dim ntb As Integer Dim ielm As Integer Dim jelm As Integer Dim kelm As Integer Dim ivx As Integer Dim ips As Integer Dim ipg As Integer Dim iva As Integer Dim ivb As Integer Dim iena(,) As Integer Dim ienb(,) As Integer Dim jeea(,) As Integer Dim jeeb(,) As Integer Dim nsra() As Integer Dim nsrb() As Integer Dim ihen(,) As Integer Dim jhen() As Integer Dim iad() As Integer Dim jstack() As Integer ReDim iena(LTE, 3) ReDim ienb(LTE, 3) ReDim jeea(LTE, 3) ReDim jeeb(LTE, 3) ReDim nsra(LTE + 2) ReDim nsrb(LTE + 2) ReDim ihen(2 * LTE + 1, 2) ReDim jhen(2 * LTE + 1) ReDim iad(2 * LTE + 1) ReDim jstack(LTE + 2) If iv = 2 Then Call EDGE(kv(1), kv(2), jac, ia) Call EDGE(kv(2), kv(1), jac, ja) ir = jac(kv(1), (ia Mod 3) + 1) jr = jac(kv(2), (ja Mod 3) + 1) mtj(kv(1), (ia Mod 3) + 1) = iq jac(kv(1), ia) = jr jac(kv(1), (ja Mod 3) + 1) = kv(2) mtj(kv(2), (ja Mod 3) + 1) = ip jac(kv(2), ja) = ir jac(kv(2), (ja Mod 3) + 1) = kv(1) If ir <> 0 Then Call EDGE(ir, kv(1), jac, iedge) jac(ir, iedge) = kv(2) End If If jr <> 0 Then Call EDGE(jr, kv(2), jac, iedge) jac(jr, iedge) = kv(1) End If iv1 = mtj(kv(1), ia) iv2 = mtj(kv(2), ja) Call DECR(iv1, kv(2), jnb, nei) Call DECR(iv2, kv(1), jnb, nei) Call INCR(iq, kv(1), jnb, nei) Call INCR(ip, kv(2), jnb, nei) Else npa = 0 npb = 0 For i = 1 To LTE + 2 nsra(i) = 0 nsrb(i) = 0 Next i For i = 1 To NELM map(i) = 0 Next i For i = 1 To iv map(kv(1)) = 1 Next i For i = 1 To iv ielm = kv(i) For j = 1 To 3 ivx = mtj(ielm, j) Call DECR(ivx, ielm, jnb, nei) Next j Next i Call PICK(iq, ip, iv, kv, mtj, jac, map, npa, npb, nsra, nsrb) Call SUBDIV(npa, nsra, px, py, nta, iena, jeea, ihen, jhen, iad, jstack) Call SUBDIV(npb, nsrb, px, py, ntb, ienb, jeeb, ihen, jhen, iad, jstack) If iv <> nta + ntb Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB POLY", "Error終了") Me.Close() End If For i = 1 To iv For j = 1 To 3 jac(kv(i), j) = 0 Next j Next i For i = 1 To nta ielm = kv(i) For j = 1 To 3 mtj(ielm, j) = iena(i, j) If jeea(i, j) <> 0 Then jac(ielm, j) = kv(jeea(i, j)) End If Next j Next i For i = 1 To ntb ielm = kv(nta + i) For j = 1 To 3 mtj(ielm, j) = ienb(i, j) If jeeb(i, j) <> 0 Then jac(ielm, j) = kv(nta + jeeb(i, j)) End If Next j Next i For i = 1 To iv ielm = kv(i) For j = 1 To 3 ivx = mtj(ielm, j) Call INCR(ivx, ielm, jnb, nei) Next j Next i For i = 1 To iv ielm = kv(i) For j = 1 To 3 jelm = jac(ielm, j) If jelm <> 0 Then GoTo poly130 ips = mtj(ielm, j) ipg = mtj(ielm, (j Mod 3) + 1) For k = 1 To jnb(ipg) kelm = nei(ipg, k) For l = 1 To 3 iva = mtj(kelm, l) ivb = mtj(kelm, (l Mod 3) + 1) If (iva = ipg) And (ivb = ips) Then jac(ielm, j) = kelm jac(kelm, l) = ielm GoTo poly130 End If Next l Next k poly130: Next j Next i End If End Sub Private Sub PICK(ByRef iq As Integer, ByRef ip As Integer, ByRef iv As Integer, ByRef kv() As Integer, _ ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef map() As Integer, _ ByRef npa As Integer, ByRef npb As Integer, ByRef nsra() As Integer, ByRef nsrb() As Integer) Dim i As Integer Dim ivx As Integer Dim jvx As Integer Dim jelm As Integer npa = 1 nsra(npa) = ip ivx = IVERT(kv(1), ip, mtj) npa = npa + 1 nsra(npa) = mtj(kv(1), (ivx Mod 3) + 1) For i = 2 To iv - 1 jvx = IVERT(kv(i), nsra(npa), mtj) jelm = jac(kv(i), jvx) If jelm = 0 Then npa = npa + 1 nsra(npa) = mtj(kv(i), (jvx Mod 3) + 1) ElseIf map(jelm) = 0 Then npa = npa + 1 nsra(npa) = mtj(kv(i), (jvx Mod 3) + 1) End If Next i npa = npa + 1 nsra(npa) = iq npb = 1 nsrb(npb) = iq ivx = IVERT(kv(iv), iq, mtj) npb = npb + 1 nsrb(npb) = mtj(kv(iv), (ivx Mod 3) + 1) For i = iv - 1 To 2 Step -1 jvx = IVERT(kv(i), nsrb(npb), mtj) jelm = jac(kv(i), jvx) If jelm = 0 Then npb = npb + 1 nsrb(npb) = mtj(kv(i), (jvx Mod 3) + 1) ElseIf map(jelm) = 0 Then npb = npb + 1 nsrb(npb) = mtj(kv(i), (jvx Mod 3) + 1) End If Next i npb = npb + 1 nsrb(npb) = ip If iv <> npa + npb - 4 Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB PICK", "Error終了") Me.Close() End If End Sub Private Sub SUBDIV(ByRef npl As Integer, ByRef nsr() As Integer, ByRef px() As Double, ByRef py() As Double, _ ByRef nte As Integer, ByRef ien(,) As Integer, ByRef jee(,) As Integer, _ ByRef ihen(,) As Integer, ByRef jhen() As Integer, ByRef iad() As Integer, ByRef jstack() As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim nnpl As Integer Dim nbs As Integer Dim ia As Integer Dim ib As Integer Dim ic As Integer Dim ix As Integer Dim xa As Double Dim ya As Double Dim xb As Double Dim yb As Double Dim see As Double nte = 0 nnpl = npl For i = 1 To LTE For j = 1 To 3 ien(i, j) = 0 jee(i, j) = 0 Next j Next i subdiv20: If 3 <= nnpl Then nbs = 1 subdiv30: If nbs <= nnpl - 1 Then ia = nsr(nbs) ib = nsr(nbs + 1) ic = nsr(((nbs + 1) Mod nnpl) + 1) xa = px(ib) - px(ia) ya = py(ib) - py(ia) xb = px(ic) - px(ia) yb = py(ic) - py(ia) see = xa * yb - xb * ya If 0.000000000000001 < see Then nte = nte + 1 ien(nte, 1) = ia ien(nte, 2) = ib ien(nte, 3) = ic nnpl = nnpl - 1 For i = nbs + 1 To nnpl nsr(i) = nsr(i + 1) Next i End If nbs = nbs + 1 GoTo subdiv30 End If GoTo subdiv20 End If ix = 0 For i = 1 To 2 * LTE + 1 ihen(i, 1) = 0 ihen(i, 2) = 0 jhen(i) = 0 iad(i) = 0 Next i For i = 1 To nte For j = 1 To 3 ia = ien(i, j) ib = ien(i, (j Mod 3) + 1) For k = 1 To ix If (ihen(k, 1) = ib) And (ihen(k, 2) = ia) Then jee(i, j) = jhen(k) jee(jhen(k), iad(k)) = i GoTo subdiv70 End If Next k ix = ix + 1 jhen(ix) = i iad(ix) = j ihen(ix, 1) = ia ihen(ix, 2) = ib subdiv70: Next j Next i Call LAWSON(nte, ien, jee, npl, px, py, jstack) End Sub Private Sub LAWSON(ByRef nte As Integer, ByRef ien(,) As Integer, ByRef jee(,) As Integer, ByRef npl As Integer, _ ByRef px() As Double, ByRef py() As Double, ByRef jstack() As Integer) Dim i As Integer Dim j As Integer Dim itop As Integer Dim maxstk As Integer Dim ncount As Integer Dim ielm As Integer Dim il As Integer Dim ir As Integer Dim jl1 As Integer Dim jl2 As Integer Dim jl3 As Integer Dim jr1 As Integer Dim jr2 As Integer Dim jr3 As Integer Dim iv1 As Integer Dim iv2 As Integer Dim iv3 As Integer Dim iv4 As Integer Dim ia As Integer Dim ib As Integer Dim iswap As Integer Dim iedge As Integer Dim xx As Double Dim yy As Double itop = 0 maxstk = npl ncount = 0 For i = 1 To nte ielm = i itop = itop + 1 jstack(itop) = IPUSH(ielm, maxstk, itop, jstack) Next i lawson20: If 0 < itop Then ncount = ncount + 1 If LTE < ncount Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB LAWSON", "Error終了") Me.Close() End If il = jstack(itop) itop = itop - 1 For j = 1 To 3 jl1 = j jl2 = (jl1 Mod 3) + 1 jl3 = (jl2 Mod 3) + 1 ir = jee(il, jl1) If ir = 0 Then GoTo lawson30 iv1 = ien(il, jl1) iv2 = ien(il, jl2) iv3 = ien(il, jl3) xx = px(ien(il, jl3)) yy = py(ien(il, jl3)) Call EDGE(ir, il, jee, jr1) jr2 = (jr1 Mod 3) + 1 jr3 = (jr2 Mod 3) + 1 iv4 = ien(ir, jr3) Call SWAP(px(iv2), py(iv2), px(iv1), py(iv1), px(iv4), py(iv4), xx, yy, iswap) If iswap = 1 Then ia = jee(il, jl2) ib = jee(ir, jr2) ien(il, jl2) = iv4 jee(il, jl1) = ib jee(il, jl2) = ir ien(ir, jr2) = iv3 jee(ir, jr1) = ia jee(ir, jr2) = il If ia <> 0 Then Call EDGE(ia, il, jee, iedge) jee(ia, iedge) = ir End If If ib <> 0 Then Call EDGE(ib, ir, jee, iedge) jee(ib, iedge) = il End If itop = itop + 1 jstack(itop) = IPUSH(il, maxstk, itop, jstack) GoTo lawson20 End If lawson30: Next j GoTo lawson20 End If End Sub Private Sub FINE(ByRef NODE As Integer, ByRef px() As Double, ByRef py() As Double, ByRef jnb() As Integer, ByRef nei(,) As Integer, _ ByRef ifix() As Integer, ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, _ ByRef idm() As Integer, ByRef istack() As Integer, ByRef angl() As Double, ByRef nedg() As Integer, _ ByRef kstack() As Integer, ByRef map() As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim l As Integer Dim ii As Integer Dim kk As Integer Dim mm As Integer Dim xx As Double Dim gx As Double Dim gy As Double Dim ar As Double Dim s As Double Dim xc As Double Dim yc As Double Dim cgrax As Double Dim cgray As Double Dim xkp As Double Dim ykp As Double Dim drate As Double Dim ntotal As Integer Dim almt As Double Dim neib As Integer Dim ielm As Integer Dim jelm As Integer Dim kelm As Integer Dim lelm As Integer Dim melm As Integer Dim j1 As Integer Dim j2 As Integer Dim j3 As Integer Dim nn As Integer Dim ls As Integer Dim le As Integer Dim lm As Integer Dim ied As Integer Dim jed As Integer Dim iv1 As Integer Dim iv2 As Integer Dim jedg As Integer ntotal = NODT If ntotal <= NODE Then ntotal = NODE almt = 1.0 / CDbl(ntotal) fine20: almt = 0.5 * almt neib = 0 For i = 1 To KTE angl(i) = 0.0 nedg(i) = 0 kstack(i) = 0 Next i For i = 1 To NELM If idm(i) <> 0 Then ielm = i neib = neib + 1 kstack(neib) = ielm Call DISTOR(ielm, mtj, px, py, drate, jedg) s = AREA(ielm, mtj, px, py) angl(ielm) = drate nedg(ielm) = jedg End If Next i For i = 1 To neib - 1 k = kstack(i) xx = angl(k) ii = i + 1 l = i For j = ii To neib kk = kstack(j) If angl(kk) <= xx Then GoTo fine60 xx = angl(kk) l = j fine60: Next j mm = kstack(i) kstack(i) = kstack(l) kstack(l) = mm Next i fine70: If neib = 0 Then GoTo fine20 If NODE < ntotal Then ielm = kstack(1) ied = nedg(ielm) jelm = jac(ielm, ied) Call EDGE(jelm, ielm, jac, jed) iv1 = mtj(ielm, ied) iv2 = mtj(ielm, (ied Mod 3) + 1) NODE = NODE + 1 px(NODE) = 0.5 * (px(iv1) + px(iv2)) py(NODE) = 0.5 * (py(iv1) + py(iv2)) Call REMESH(ielm, ied, jelm, jed, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, istack) If idm(ielm) <> idm(jelm) Then ifix(NODE) = 1 Else gx = 0.0 gy = 0.0 ar = 0.0 For j = 1 To jnb(NODE) kelm = nei(NODE, j) j1 = mtj(kelm, 1) j2 = mtj(kelm, 2) j3 = mtj(kelm, 3) s = AREA(kelm, mtj, px, py) xc = (px(j1) + px(j2) + px(j3)) / 3.0 yc = (py(j1) + py(j2) + py(j3)) / 3.0 ar = ar + s gx = gx + s * xc gy = gy + s * yc Next j cgrax = gx / ar cgray = gy / ar xkp = px(NODE) ykp = py(NODE) px(NODE) = cgrax py(NODE) = cgray For j = 1 To jnb(NODE) kelm = nei(NODE, j) s = AREA(kelm, mtj, px, py) If s <= 0.1 * almt Then px(NODE) = xkp py(NODE) = ykp Exit For End If Next j End If For i = 1 To NELM map(i) = 1 Next i For i = 1 To jnb(NODE) lelm = nei(NODE, i) map(lelm) = 0 angl(lelm) = 0.0 nedg(lelm) = 0 Next i nn = 0 For i = 1 To neib If map(kstack(i)) = 1 Then nn = nn + 1 kstack(nn) = kstack(i) End If Next i neib = nn For i = 1 To jnb(NODE) melm = nei(NODE, i) s = AREA(melm, mtj, px, py) If (almt < s) And (idm(melm) <> 0) Then Call DISTOR(melm, mtj, px, py, drate, jedg) angl(melm) = drate nedg(melm) = jedg If angl(kstack(1)) < drate Then For j = 1 To neib kstack(neib - j + 2) = kstack(neib - j + 1) Next j neib = neib + 1 kstack(1) = melm ElseIf drate < angl(kstack(neib)) Then neib = neib + 1 kstack(neib) = melm Else ls = 1 le = neib fine160: If le - ls <> 1 Then lm = CInt((ls + le) / 2) If angl(kstack(lm)) < drate Then le = lm Else ls = lm End If GoTo fine160 End If For k = neib To le Step -1 kstack(k + 1) = kstack(k) Next k neib = neib + 1 kstack(le) = melm End If End If Next i GoTo fine70 End If End Sub Private Sub REMESH(ByRef ielm As Integer, ByRef ied As Integer, ByRef jelm As Integer, ByRef jed As Integer, ByRef NODE As Integer, _ ByRef px() As Double, ByRef py() As Double, ByRef jnb() As Integer, ByRef nei(,) As Integer, _ ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef idm() As Integer, ByRef istack() As Integer) Dim itop As Integer Dim maxstk As Integer Dim ip As Integer Dim xp As Double Dim yp As Double Dim ibj As Integer Dim jbj As Integer Dim ia As Integer Dim ib As Integer Dim ic As Integer Dim id As Integer Dim iv1 As Integer Dim iv2 As Integer Dim iv3 As Integer Dim iv4 As Integer Dim iedge As Integer Dim il As Integer Dim ir As Integer Dim iera As Integer Dim ierb As Integer Dim ierl As Integer Dim iswap As Integer itop = 0 maxstk = NODE ip = NODE xp = px(ip) yp = py(ip) ibj = (ied Mod 3) + 1 jbj = (jed Mod 3) + 1 ia = jac(ielm, ibj) ib = jac(ielm, (ibj Mod 3) + 1) ic = jac(jelm, jbj) id = jac(jelm, (jbj Mod 3) + 1) iv1 = mtj(ielm, ibj) iv2 = mtj(ielm, (ibj Mod 3) + 1) iv3 = mtj(jelm, jbj) iv4 = mtj(jelm, (jbj Mod 3) + 1) mtj(ielm, 1) = ip mtj(ielm, 2) = iv1 mtj(ielm, 3) = iv2 jac(ielm, 1) = NELM + 2 jac(ielm, 2) = ia jac(ielm, 3) = NELM + 1 NELM = NELM + 1 idm(NELM) = idm(ielm) mtj(NELM, 1) = ip mtj(NELM, 2) = iv2 mtj(NELM, 3) = iv3 jac(NELM, 1) = ielm jac(NELM, 2) = ib jac(NELM, 3) = jelm mtj(jelm, 1) = ip mtj(jelm, 2) = iv3 mtj(jelm, 3) = iv4 jac(jelm, 1) = NELM jac(jelm, 2) = ic jac(jelm, 3) = NELM + 1 NELM = NELM + 1 idm(NELM) = idm(jelm) mtj(NELM, 1) = ip mtj(NELM, 2) = iv4 mtj(NELM, 3) = iv1 jac(NELM, 1) = jelm jac(NELM, 2) = id jac(NELM, 3) = ielm If iv1 <= KTJ Then nei(iv1, NEIBOR(iv1, jelm, jnb, nei)) = NELM End If Call INCR(iv2, NELM - 1, jnb, nei) If iv3 <= KTJ Then nei(iv3, NEIBOR(iv3, ielm, jnb, nei)) = NELM - 1 End If Call INCR(iv4, NELM, jnb, nei) jnb(ip) = 4 nei(ip, 1) = ielm nei(ip, 2) = NELM - 1 nei(ip, 3) = jelm nei(ip, 4) = NELM If ia <> 0 Then If idm(ia) = idm(ielm) Then itop = itop + 1 istack(itop) = IPUSH(ielm, maxstk, itop, istack) End If End If If ib <> 0 Then Call EDGE(ib, ielm, jac, iedge) jac(ib, iedge) = NELM - 1 If idm(ib) = idm(NELM - 1) Then itop = itop + 1 istack(itop) = IPUSH(NELM - 1, maxstk, itop, istack) End If End If If ic <> 0 Then If idm(ic) = idm(jelm) Then itop = itop + 1 istack(itop) = IPUSH(jelm, maxstk, itop, istack) End If End If If id <> 0 Then Call EDGE(id, jelm, jac, iedge) jac(id, iedge) = NELM If idm(id) = idm(NELM) Then itop = itop + 1 istack(itop) = IPUSH(NELM, maxstk, itop, istack) End If End If remesh10: If 0 < itop Then il = istack(itop) itop = itop - 1 ir = jac(il, 2) Call EDGE(ir, il, jac, ierl) iera = (ierl Mod 3) + 1 ierb = (iera Mod 3) + 1 iv1 = mtj(ir, ierl) iv2 = mtj(ir, iera) iv3 = mtj(ir, ierb) Call SWAP(px(iv1), py(iv1), px(iv2), py(iv2), px(iv3), py(iv3), xp, yp, iswap) If iswap = 1 Then ia = jac(ir, iera) ib = jac(ir, ierb) ic = jac(il, 3) mtj(il, 3) = iv3 jac(il, 2) = ia jac(il, 3) = ir mtj(ir, 1) = ip mtj(ir, 2) = iv3 mtj(ir, 3) = iv1 jac(ir, 1) = il jac(ir, 2) = ib jac(ir, 3) = ic Call DECR(iv1, il, jnb, nei) Call DECR(iv2, ir, jnb, nei) Call INCR(ip, ir, jnb, nei) Call INCR(iv3, il, jnb, nei) If ia <> 0 Then Call EDGE(ia, ir, jac, iedge) jac(ia, iedge) = il If idm(ia) = idm(il) Then itop = itop + 1 istack(itop) = IPUSH(il, maxstk, itop, istack) End If End If If ib <> 0 Then If idm(ib) = idm(ir) Then itop = itop + 1 istack(itop) = IPUSH(ir, maxstk, itop, istack) End If End If If ic <> 0 Then Call EDGE(ic, il, jac, iedge) jac(ic, iedge) = ir End If End If GoTo remesh10 End If End Sub Private Sub LAPLAS(ByRef NODE As Integer, ByRef ifix() As Integer, ByRef mtj(,) As Integer, ByRef px() As Double, ByRef py() As Double, _ ByRef jnb() As Integer, ByRef nei(,) As Integer) Dim itera As Integer = 5 Dim it As Integer Dim i As Integer Dim j As Integer Dim ielm As Integer Dim j1 As Integer Dim j2 As Integer Dim j3 As Integer Dim gx As Double Dim gy As Double Dim ar As Double Dim s As Double Dim xc As Double Dim yc As Double Dim cgrax As Double Dim cgray As Double For it = 1 To itera For i = 1 To NODE If ifix(i) = 0 Then gx = 0.0 gy = 0.0 ar = 0.0 For j = 1 To jnb(i) ielm = nei(i, j) j1 = mtj(ielm, 1) j2 = mtj(ielm, 2) j3 = mtj(ielm, 3) s = AREA(ielm, mtj, px, py) xc = (px(j1) + px(j2) + px(j3)) / 3.0 yc = (py(j1) + py(j2) + py(j3)) / 3.0 ar = ar + s gx = gx + s * xc gy = gy + s * yc Next j cgrax = gx / ar cgray = gy / ar px(i) = cgrax py(i) = cgray End If Next i Next it End Sub Private Sub REMOVE(ByRef NEX As Integer, ByRef index() As Integer, ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef idm() As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim l As Integer Dim iz As Integer Dim inelm As Integer Dim ielm As Integer Dim jelm As Integer Dim kelm As Integer Dim lelm As Integer Dim mkp(3) As Integer Dim jkp(3) As Integer Dim ikp As Integer Dim kedg As Integer Dim ledg As Integer ielm = 0 index(1) = 1 For i = 1 To NEX iz = i inelm = 0 For j = index(i) To NELM If idm(j) = iz Then ielm = ielm + 1 jelm = j inelm = inelm + 1 If ielm <> jelm Then For k = 1 To 3 mkp(k) = mtj(ielm, k) jkp(k) = jac(ielm, k) Next k ikp = idm(ielm) For k = 1 To 3 kelm = jac(jelm, k) If kelm <> 0 Then Call EDGE(kelm, jelm, jac, kedg) jac(kelm, kedg) = ielm + KTE + 1 End If Next k For l = 1 To 3 lelm = jkp(l) If lelm <> 0 Then Call EDGE(lelm, ielm, jac, ledg) jac(lelm, ledg) = jelm + KTE + 1 End If Next l For k = 1 To 3 jac(ielm, k) = jac(ielm, k) Mod (KTE + 1) jac(jelm, k) = jac(jelm, k) Mod (KTE + 1) kelm = jac(ielm, k) lelm = jac(jelm, k) For l = 1 To 3 jac(kelm, l) = jac(kelm, l) Mod (KTE + 1) jac(lelm, l) = jac(lelm, l) Mod (KTE + 1) Next l Next k For k = 1 To 3 jkp(k) = jac(ielm, k) Next k For k = 1 To 3 mtj(ielm, k) = mtj(jelm, k) jac(ielm, k) = jac(jelm, k) mtj(jelm, k) = mkp(k) jac(jelm, k) = jkp(k) Next k idm(ielm) = idm(jelm) idm(jelm) = ikp End If End If Next j index(i + 1) = index(i) + inelm Next i For i = 1 To ielm For j = 1 To 3 If ielm < jac(i, j) Then jac(i, j) = 0 End If Next j Next i For i = ielm + 1 To NELM For j = 1 To 3 mtj(i, j) = 0 jac(i, j) = 0 Next j Next i NELM = ielm End Sub Private Sub CHECK(ByRef NELM As Integer, ByRef mtj(,) As Integer, ByRef jac(,) As Integer) Dim i As Integer Dim j As Integer Dim k As Integer Dim ielm As Integer Dim jelm As Integer Dim kelm As Integer Dim ia As Integer Dim ib As Integer Dim ja As Integer Dim jb As Integer For i = 1 To NELM For j = 1 To 3 ielm = i ia = mtj(i, j) ib = mtj(i, (j Mod 3) + 1) jelm = jac(i, j) If jelm = 0 Then GoTo check10 For k = 1 To 3 kelm = jac(jelm, k) If kelm = ielm Then ja = mtj(jelm, k) jb = mtj(jelm, (k Mod 3) + 1) If (ia = jb) And (ib = ja) Then GoTo check10 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB CHECK", "Error終了") Me.Close() End If Next k 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB CHECK", "Error終了") Me.Close() check10: Next j Next i End Sub Private Sub DATA(ByRef NEX As Integer, ByRef NODE As Integer, ByRef NELM As Integer, ByRef mtj(,) As Integer, _ ByRef jac(,) As Integer, ByRef idm() As Integer, ByRef px() As Double, ByRef py() As Double) Dim sw As System.IO.StreamWriter Dim dat As String Dim fname2 As String Dim i As Integer fname2 = My.Forms.Form1.Label2.Text sw = New System.IO.StreamWriter(fname2, False, System.Text.Encoding.Default) dat = NODE.ToString & "," & NELM.ToString : sw.WriteLine(dat) For i = 1 To NELM dat = mtj(i, 1).ToString & "," & mtj(i, 2).ToString & "," & mtj(i, 3).ToString dat = dat & "," & idm(i).ToString sw.WriteLine(dat) Next i For i = 1 To NODE dat = px(i).ToString("0.000") & "," & py(i).ToString("0.000") sw.WriteLine(dat) Next i sw.Close() End Sub Private Sub INCR(ByRef n As Integer, ByRef l As Integer, ByRef jnb() As Integer, ByRef nei(,) As Integer) If n <= KTJ Then jnb(n) = jnb(n) + 1 If KCM < jnb(n) Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB INCR", "Error終了") Me.Close() End If nei(n, jnb(n)) = l End If End Sub Private Sub DECR(ByRef n As Integer, ByRef l As Integer, ByRef jnb() As Integer, ByRef nei(,) As Integer) Dim i As Integer Dim inb As Integer If n <= KTJ Then inb = NEIBOR(n, l, jnb, nei) For i = inb To jnb(n) - 1 nei(n, i) = nei(n, i + 1) Next i nei(n, jnb(n)) = 0 jnb(n) = jnb(n) - 1 If jnb(n) = 0 Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN SUB DECR", "Error終了") Me.Close() End If End If End Sub Private Sub DISTOR(ByRef n As Integer, ByRef mtj(,) As Integer, ByRef px() As Double, ByRef py() As Double, ByRef drate As Double, ByRef jedg As Integer) Dim dst(3) As Double Dim xa As Double Dim ya As Double Dim xb As Double Dim yb As Double Dim xc As Double Dim yc As Double Dim ang1 As Double Dim ang2 As Double Dim ang3 As Double xa = px(mtj(n, 1)) ya = py(mtj(n, 1)) xb = px(mtj(n, 2)) yb = py(mtj(n, 2)) xc = px(mtj(n, 3)) yc = py(mtj(n, 3)) ang1 = THETA(xc, yc, xa, ya, xb, yb) ang2 = THETA(xa, ya, xb, yb, xc, yc) ang3 = pi - ang1 - ang2 dst(1) = ang1 - pi / 3.0 dst(2) = ang2 - pi / 3.0 dst(3) = ang3 - pi / 3.0 drate = Math.Abs(dst(1)) + Math.Abs(dst(2)) + Math.Abs(dst(3)) jedg = 1 If dst(jedg) < dst(2) Then jedg = 2 If dst(jedg) < dst(3) Then jedg = 3 End Sub Private Sub LOCATE(ByRef xp As Double, ByRef yp As Double, ByRef px() As Double, ByRef py() As Double, _ ByRef mtj(,) As Integer, ByRef jac(,) As Integer, ByRef NELM As Integer, ByRef itri As Integer) Dim i As Integer Dim ia As Integer Dim ib As Integer itri = NELM locate10: For i = 1 To 3 ia = mtj(itri, i) ib = mtj(itri, (i Mod 3) + 1) If (px(ia) - xp) * (py(ib) - yp) < (py(ia) - yp) * (px(ib) - xp) Then itri = jac(itri, i) GoTo locate10 End If Next i End Sub Private Sub SWAP(ByRef x1 As Double, ByRef y1 As Double, ByRef x2 As Double, ByRef y2 As Double, ByRef x3 As Double, ByRef y3 As Double, _ ByRef xp As Double, ByRef yp As Double, ByRef iswap As Integer) Dim x13 As Double Dim y13 As Double Dim x23 As Double Dim y23 As Double Dim x1p As Double Dim y1p As Double Dim x2p As Double Dim y2p As Double Dim cosa As Double Dim cosb As Double Dim sina As Double Dim sinb As Double x13 = x1 - x3 y13 = y1 - y3 x23 = x2 - x3 y23 = y2 - y3 x1p = x1 - xp y1p = y1 - yp x2p = x2 - xp y2p = y2 - yp cosa = x13 * x23 + y13 * y23 cosb = x2p * x1p + y1p * y2p If (0.0 <= cosa) And (0.0 <= cosb) Then iswap = 0 ElseIf (cosa < 0.0) And (cosb < 0.0) Then iswap = 1 Else sina = x13 * y23 - x23 * y13 sinb = x2p * y1p - x1p * y2p If (sina * cosb + sinb * cosa) < 0 Then iswap = 1 Else iswap = 0 End If End If End Sub Private Sub EDGE(ByRef l As Integer, ByRef k As Integer, ByRef jac(,) As Integer, ByRef iedge As Integer) Dim i As Integer iedge = -1 For i = 1 To 3 If jac(l, i) = k Then iedge = i Exit For End If Next i 'ERRORメッセージ&実行中止 If iedge = -1 Then MessageBox.Show("ERROR IN SUB EDGE", "Error終了") Me.Close() End If End Sub Private Function IPUSH(ByRef item As Integer, ByRef maxstk As Integer, ByRef itop As Integer, ByRef istack() As Integer) As Integer If maxstk < itop Then 'ERRORメッセージ&実行中止 MessageBox.Show("ERROR IN FUNCTION IPUSH", "Error終了") Me.Close() Else IPUSH = item End If End Function Private Function NEIBOR(ByRef n As Integer, ByRef l As Integer, ByRef jnb() As Integer, ByRef nei(,) As Integer) As Integer Dim i As Integer NEIBOR = -1 For i = 1 To jnb(n) If nei(n, i) = l Then NEIBOR = i Exit For End If Next i 'ERRORメッセージ&実行中止 If NEIBOR = -1 Then MessageBox.Show("ERROR IN FUNCTION NEIBOR", "Error終了") Me.Close() End If End Function Private Function IVERT(ByRef l As Integer, ByRef k As Integer, ByRef mtj(,) As Integer) As Integer Dim i As Integer IVERT = -1 For i = 1 To 3 If mtj(l, i) = k Then IVERT = i Exit For End If Next i 'ERRORメッセージ&実行中止 If IVERT = -1 Then MessageBox.Show("ERROR IN FUNCTION IVERT", "Error終了") Me.Close() End If End Function Private Function AREA(ByRef ielm As Integer, ByRef mtj(,) As Integer, ByRef px() As Double, ByRef py() As Double) As Double Dim j1 As Integer Dim j2 As Integer Dim j3 As Integer j1 = mtj(ielm, 1) j2 = mtj(ielm, 2) j3 = mtj(ielm, 3) AREA = 0.5 * (px(j1) * py(j2) + px(j2) * py(j3) + px(j3) * py(j1) - px(j1) * py(j3) - px(j2) * py(j1) - px(j3) * py(j2)) End Function Private Function THETA(ByRef x0 As Double, ByRef y0 As Double, ByRef x1 As Double, ByVal y1 As Double, _ ByRef x2 As Double, ByRef y2 As Double) As Double Dim err As Double Dim xa As Double Dim ya As Double Dim xb As Double Dim yb As Double Dim prdin As Double Dim prdex As Double err = 0.000000000000001 xa = x1 - x0 ya = y1 - y0 xb = x2 - x0 yb = y2 - y0 prdin = xa * xb + ya * yb prdex = xa * yb - xb * ya If Math.Abs(prdin) < err Then THETA = pi / 2.0 Else THETA = Math.Atan(prdex / prdin) If THETA < 0.0 Then THETA = THETA + pi End If If pi < THETA Then THETA = THETA - pi End If End If End Function End Class