module defparam implicit none integer,parameter::KBD = 500 !最大境界数 integer,parameter::KTJ = 10000 !最大節点数 integer,parameter::KCM = 100 !節点に集まる最大要素数 integer,parameter::LTE = 100 !多角形の最大辺数 integer,parameter::KTE = 2 * KTJ + 1 integer::INPNODT !指定節点数(入力データ) real(4),parameter::pi=3.14159265358979323846 end module defparam !----------------------------------------------------------------------------------------- program f90_FEM_MESH3 !----------------------------------------------------------------------------------------- use defparam implicit none integer::ibex(1:KBD) =0 !個々の外部境界上の節点数 integer::ibin(1:KBD) =0 !個々の内部境界上の節点数 integer::ibno(1:KBD,1:KTJ) =0 !外部境界上の節点番号 integer::mindex(1:KBD+1) =0 !作業領域 real(4)::px(1:KTJ+3) =0.0 !節点のx座標 real(4)::py(1:KTJ+3) =0.0 !節点のy座標 integer::mtj(1:2*KTJ+1,1:3) =0 !要素-節点関係 integer::jac(1:2*KTJ+1,1:3) =0 !要素-隣接要素関係 integer::idm(1:2*KTJ+1) =0 !要素-領域関係 integer::NEX =0 !外部境界数 integer::NIN =0 !内部境界数 integer::NOB =0 !境界上の節点総数 integer::NIB =0 !境界内部の節点総数 integer::NODE =0 !節点総数 integer::NELM =0 !要素総数 character(len=50)::fnameR !入力ファイル名 character(len=50)::fnameW !出力ファイル名 call getarg(1,fnameR) call getarg(2,fnameW) call INPUT(NEX, NIN, ibex, ibin, ibno, NOB, NIB, px, py, fnameR) call MODEL(NEX, NIN, ibex, ibin, ibno, NOB, NIB, mindex, NODE, px, py, NELM, mtj, jac, idm) call DATA(NODE, NELM, mtj, jac, idm, px, py, fnameW) stop contains !----------------------------------------------------------------------------------------- subroutine INPUT(NEX, NIN, ibex, ibin, ibno, NOB, NIB, px, py, fnameR) !----------------------------------------------------------------------------------------- integer,intent(inout)::NEX integer,intent(inout)::NIN integer,intent(inout)::ibex(1:KBD) integer,intent(inout)::ibin(1:KBD) integer,intent(inout)::ibno(1:KBD,1:KTJ) integer,intent(inout)::NOB integer,intent(inout)::NIB real(4),intent(inout)::px(1:KTJ+3) real(4),intent(inout)::py(1:KTJ+3) character(len=50)::fnameR integer::i,j open(11,file=fnameR,status='old') read(11,*) INPNODT read(11,*) NEX,NIN if(0 < NEX)then do i = 1,NEX read(11,*) ibex(i) end do end if if(0 < NIN)then do i = 1,NIN read(11,*) ibin(i) end do end if do i = 1,NEX read(11,*) (ibno(i,j),j=1,ibex(i)) end do read(11,*) NOB,NIB do i = 1,NOB + NIB read(11,*) px(i),py(i) end do close(11) end subroutine INPUT !----------------------------------------------------------------------------------------- subroutine MODEL(NEX, NIN, ibex, ibin, ibno, NOB, NIB, mindex, NODE, px, py, NELM, mtj, jac, idm) !----------------------------------------------------------------------------------------- integer,intent(inout)::NEX integer,intent(inout)::NIN integer,intent(inout)::ibex(1:KBD) integer,intent(inout)::ibin(1:KBD) integer,intent(inout)::ibno(1:KBD,1:KTJ) integer,intent(inout)::NOB integer,intent(inout)::NIB integer,intent(inout)::mindex(1:KBD+1) integer,intent(inout)::NODE real(4),intent(inout)::px(1:KTJ+3) real(4),intent(inout)::py(1:KTJ+3) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer::i,ntp,ia,ib,ic real(4)::xmin,xmax,ymin,ymax,rax,ray,dmax integer::jnb(1:KTJ) =0 integer::nei(1:KTJ,1:KCM) =0 integer::ifix(1:KTJ) =0 real(4)::angl(1:2*KTJ+1) =0.0 integer::nedg(1:2*KTJ+1) =0 integer::kstack(1:2*KTJ+1) =0 integer::list(1:KTJ) =0 integer::istack(1:KTJ) =0 integer::kv(1:2*KTJ+1) =0 integer::ibr(1:2*KTJ+1) =0 integer::iadres(1:KTJ+3) =0 integer::map(1:2*KTJ+1) =0 ntp = NOB + NIB xmin = px(1) xmax = xmin ymin = py(1) ymax = ymin do i = 2,ntp if(px(i) < xmin)xmin = px(i) if(xmax < px(i))xmax = px(i) if(py(i) < ymin)ymin = py(i) if(ymax < py(i))ymax = py(i) end do rax = xmax - xmin ray = ymax - ymin dmax = rax if(dmax < ray)dmax = ray do i = 1,ntp px(i) = (px(i) - xmin) / dmax py(i) = (py(i) - ymin) / dmax end do 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) do i = 1,NODE ifix(i) = 1 end do call FINE(NODE, px, py, jnb, nei, ifix, NELM, mtj, jac, idm, istack, angl, nedg, kstack, map) if(NELM /= 2 * NODE + 1)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB MODEL: stop' stop end if call LAPLAS(NODE, ifix, mtj, px, py, jnb, nei) call REMOVE(NEX, mindex, NELM, mtj, jac, idm) call CHECK(NELM, mtj, jac) do i = 1,NODE px(i) = px(i) * dmax + xmin py(i) = py(i) * dmax + ymin end do end subroutine MODEL !----------------------------------------------------------------------------------------- subroutine ROUGH(NEX, NIN, ibex, ibin, ibno, ntp, NIB, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, & & list, iadres, istack, kv, ibr, map) !----------------------------------------------------------------------------------------- integer,intent(in)::NEX integer,intent(in)::NIN integer,intent(in)::ibex(1:KBD) integer,intent(in)::ibin(1:KBD) integer,intent(in)::ibno(1:KBD,1:KTJ) integer,intent(in)::ntp integer,intent(in)::NIB integer,intent(inout)::NODE real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer,intent(inout)::list(1:KTJ) integer,intent(inout)::iadres(1:KTJ+3) integer,intent(inout)::istack(1:KTJ) integer,intent(inout)::kv(1:2*KTJ+1) integer,intent(inout)::ibr(1:2*KTJ+1) integer,intent(inout)::map(1:2*KTJ+1) integer::i,j,k,nb,np,ma,mb,mc,ms,mp,js,jz,loc,it real(4)::xs,ys,xa,ya do i = 1,NEX jz = 0 nb = i np = ibex(nb) do j = 1,np list(j) = ibno(nb, j) end do call BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE) do k = 1,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 end do end do do i = 1,NIN nb = i np = ibin(nb) do j = 1,np list(j) = NODE + j end do 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) do k = 1,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 end do end do do i = 1,KTJ + 3 iadres(i) = 0 end do do i = 1,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) end do if(NODE /= ntp)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB ROUGH: stop' stop end if end subroutine ROUGH !----------------------------------------------------------------------------------------- subroutine BOUGEN(jz, np, list, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack, kv, ibr, map, NODE) !----------------------------------------------------------------------------------------- integer,intent(in)::jz integer,intent(in)::np integer,intent(in)::list(1:KTJ) integer,intent(in)::ntp real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer,intent(inout)::iadres(1:KTJ+3) integer,intent(inout)::istack(1:KTJ) integer,intent(inout)::kv(1:2*KTJ+1) integer,intent(inout)::ibr(1:2*KTJ+1) integer,intent(inout)::map(1:2*KTJ+1) integer,intent(inout)::NODE integer::i,j,k,inp,js,ip,iq,iv inp = 0 do i = 1,KTJ + 3 iadres(i) = 0 end do do i = 1,np iadres(list(i)) = i end do 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 loopi: do i = 1,np ip = list(mod(i,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 do j = 1,jnb(iq) do k = 1,3 if(mtj(nei(iq, j), k) == ip)cycle loopi end do end do 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) end do loopi NODE = NODE + inp end subroutine BOUGEN !----------------------------------------------------------------------------------------- subroutine DELAUN(jz, js, jg, ntp, px, py, jnb, nei, NELM, mtj, jac, idm, iadres, istack) !----------------------------------------------------------------------------------------- integer,intent(in)::jz integer,intent(in)::js integer,intent(in)::jg integer,intent(in)::ntp real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer,intent(in)::iadres(1:KTJ+3) integer,intent(inout)::istack(1:KTJ) integer::i,itop,maxstk,ip,it,ia,ib,ic,iv1,iv2,iv3,idf,ms,iedge,il,ir,iera,ierb,ierl,iswap real(4)::xp,yp itop = 0 maxstk = ntp do i = js,jg ip = i xp = px(ip) yp = py(ip) call LOCATE(xp, yp, px, py, mtj, jac, NELM, it) if(idm(it) /= jz)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB DELAUN: stop' stop 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 = 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) 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 = 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) end if end if if(ic /= 0)then call EDGE(ic, it, jac, iedge) jac(ic, iedge) = NELM ms = iadres(iv3) * iadres(iv1) idf = 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) end if end if do if(0 < itop)then il = istack(itop) itop = itop - 1 ir = jac(il, 2) call EDGE(ir, il, jac, ierl) iera = mod(ierl,3) + 1 ierb = mod(iera,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 = 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) end if end if if(ib /= 0)then ms = iadres(iv3) * iadres(iv1) idf = 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) end if end if if(ic /= 0)then call EDGE(ic, il, jac, iedge) jac(ic, iedge) = ir end if end if else exit end if end do end do end subroutine DELAUN !----------------------------------------------------------------------------------------- subroutine SEARCH(iz, ip, iq, jnb, nei, NELM, mtj, jac, idm, istack, iv, kv, iadres, ibr) !----------------------------------------------------------------------------------------- integer,intent(in)::iz integer,intent(in)::ip integer,intent(in)::iq integer,intent(in)::jnb(1:KTJ) integer,intent(in)::nei(1:KTJ,1:KCM) integer,intent(in)::NELM integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer,intent(in)::jac(1:2*KTJ+1,1:3) integer,intent(in)::idm(1:2*KTJ+1) integer,intent(inout)::istack(1:KTJ) integer,intent(inout)::iv integer,intent(inout)::kv(1:2*KTJ+1) integer,intent(in)::iadres(1:KTJ+3) integer,intent(inout)::ibr(1:2*KTJ+1) integer::i,j,ia,ib,ja,jb,mstk,nbr,idf,ms,ielm,jelm,kelm,mmin,jr iv = 0 mstk = 0 nbr = 0 do i = 1,NELM kv(i) = 0 ibr(i) = 0 end do do i = 1,KTJ istack(i) = 0 end do do i = 1,jnb(iq) nbr = nbr + 1 ibr(nei(iq, i)) = nbr end do do i = 1,jnb(iq) ielm = nei(iq, i) j = IVERT(ielm, iq, mtj) ja = mod(j,3) + 1 jb = mod(ja,3) + 1 ia = mtj(ielm, ja) ib = mtj(ielm, jb) idf = abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) if((idf == 1).and.(ms /= 0))cycle jelm = jac(ielm, ja) if(jelm == 0)cycle if(idm(jelm) /= iz)cycle nbr = nbr + 1 ibr(jelm) = nbr if(mtj(jelm, 1) == ip)goto 80 if(mtj(jelm, 2) == ip)goto 80 if(mtj(jelm, 3) == ip)goto 80 mstk = mstk + 1 istack(mstk) = jelm end do do if(mstk == 0)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB SEARCH: stop' stop end if ielm = istack(1) mstk = mstk - 1 do i = 1,mstk istack(i) = istack(i + 1) end do istack(mstk + 1) = 0 do j = 1,3 ja = mod(j,3) + 1 ia = mtj(ielm, j) ib = mtj(ielm, ja) idf = abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) if((idf == 1).and.(ms /= 0))cycle jelm = jac(ielm, j) if(jelm == 0)cycle if(ibr(jelm) /= 0)cycle if(idm(jelm) /= iz)cycle nbr = nbr + 1 ibr(jelm) = nbr if(mtj(jelm, 1) == ip)goto 80 if(mtj(jelm, 2) == ip)goto 80 if(mtj(jelm, 3) == ip)goto 80 mstk = mstk + 1 istack(mstk) = jelm end do end do 80 continue do iv = iv + 1 kv(iv) = jelm if(ibr(jelm) <= jnb(iq))exit mmin = KTE + 1 do j = 1,3 jr = jac(jelm, j) if(jr == 0)cycle if(ibr(jr) == 0)cycle ja = mod(j,3) + 1 ia = mtj(jelm, j) ib = mtj(jelm, ja) idf = abs(iadres(ia) - iadres(ib)) ms = iadres(ia) * iadres(ib) if((idf == 1).and.(ms /= 0))cycle if(ibr(jr) < mmin)then kelm = jr mmin = ibr(jr) end if end do if(mmin == KTE + 1)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB SEARCH: stop' stop end if jelm = kelm end do end subroutine SEARCH !----------------------------------------------------------------------------------------- subroutine POLY(iq, ip, iv, kv, px, py, NELM, mtj, jac, jnb, nei, map) !----------------------------------------------------------------------------------------- integer,intent(in)::iq integer,intent(in)::ip integer,intent(in)::iv integer,intent(in)::kv(1:2*KTJ+1) real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(in)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::map(1:2*KTJ+1) integer::i,j,k,l,ia,ja,ir,jr,iv1,iv2,iedge,npa,npb,nta,ntb integer::ielm,jelm,kelm,ivx,ips,ipg,iva,ivb integer::iena(1:LTE,1:3) =0 integer::ienb(1:LTE,1:3) =0 integer::jeea(1:LTE,1:3) =0 integer::jeeb(1:LTE,1:3) =0 integer::nsra(1:LTE+2) =0 integer::nsrb(1:LTE+2) =0 integer::ihen(1:2*LTE+1,1:2) =0 integer::jhen(1:2*LTE+1) =0 integer::iad(1:2*LTE+1) =0 integer::jstack(1:LTE+2) =0 if(iv == 2)then call EDGE(kv(1), kv(2), jac, ia) call EDGE(kv(2), kv(1), jac, ja) ir = jac(kv(1), mod(ia,3) + 1) jr = jac(kv(2), mod(ja,3) + 1) mtj(kv(1), mod(ia,3) + 1) = iq jac(kv(1), ia) = jr jac(kv(1), mod(ja,3) + 1) = kv(2) mtj(kv(2), mod(ja,3) + 1) = ip jac(kv(2), ja) = ir jac(kv(2), mod(ja,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 do i = 1,LTE + 2 nsra(i) = 0 nsrb(i) = 0 end do do i = 1,NELM map(i) = 0 end do do i = 1,iv map(kv(1)) = 1 end do do i = 1,iv ielm = kv(i) do j = 1,3 ivx = mtj(ielm, j) call DECR(ivx, ielm, jnb, nei) end do end do 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 !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB POLY: stop' stop end if do i = 1,iv do j = 1,3 jac(kv(i), j) = 0 end do end do do i = 1,nta ielm = kv(i) do j = 1,3 mtj(ielm, j) = iena(i, j) if(jeea(i, j) /= 0)then jac(ielm, j) = kv(jeea(i, j)) end if end do end do do i = 1,ntb ielm = kv(nta + i) do j = 1,3 mtj(ielm, j) = ienb(i, j) if(jeeb(i, j) /= 0)then jac(ielm, j) = kv(nta + jeeb(i, j)) end if end do end do do i = 1,iv ielm = kv(i) do j = 1,3 ivx = mtj(ielm, j) call INCR(ivx, ielm, jnb, nei) end do end do do i = 1,iv ielm = kv(i) loopj: do j = 1,3 jelm = jac(ielm, j) if(jelm /= 0)cycle ips = mtj(ielm, j) ipg = mtj(ielm, mod(j,3) + 1) do k = 1,jnb(ipg) kelm = nei(ipg, k) do l = 1,3 iva = mtj(kelm, l) ivb = mtj(kelm, mod(l,3) + 1) if((iva == ipg).and.(ivb == ips))then jac(ielm, j) = kelm jac(kelm, l) = ielm cycle loopj end if end do end do end do loopj end do end if end subroutine POLY !----------------------------------------------------------------------------------------- subroutine PICK(iq, ip, iv, kv, mtj, jac, map, npa, npb, nsra, nsrb) !----------------------------------------------------------------------------------------- integer,intent(in)::iq integer,intent(in)::ip integer,intent(in)::iv integer,intent(in)::kv(1:2*KTJ+1) integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer,intent(in)::jac(1:2*KTJ+1,1:3) integer,intent(in)::map(1:2*KTJ+1) integer,intent(inout)::npa integer,intent(inout)::npb integer,intent(inout)::nsra(1:LTE+2) integer,intent(inout)::nsrb(1:LTE+2) integer::i,ivx,jvx,jelm npa = 1 nsra(npa) = ip ivx = IVERT(kv(1), ip, mtj) npa = npa + 1 nsra(npa) = mtj(kv(1), mod(ivx,3) + 1) do i = 2,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), mod(jvx,3) + 1) else if(map(jelm) == 0)then npa = npa + 1 nsra(npa) = mtj(kv(i), mod(jvx,3) + 1) end if end do 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), mod(ivx,3) + 1) do i = iv - 1,2,-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), mod(jvx,3) + 1) else if(map(jelm) == 0)then npb = npb + 1 nsrb(npb) = mtj(kv(i), mod(jvx,3) + 1) end if end do npb = npb + 1 nsrb(npb) = ip if(iv /= npa + npb - 4)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB PICK: stop' stop end if end subroutine PICK !----------------------------------------------------------------------------------------- subroutine SUBDIV(npl, nsr, px, py, nte, ien, jee, ihen, jhen, iad, jstack) !----------------------------------------------------------------------------------------- integer,intent(in)::npl integer,intent(inout)::nsr(1:LTE+2) real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::nte integer,intent(inout)::ien(1:LTE,1:3) integer,intent(inout)::jee(1:LTE,1:3) integer,intent(inout)::ihen(1:2*LTE+1,1:2) integer,intent(inout)::jhen(1:2*LTE+1) integer,intent(inout)::iad(1:2*LTE+1) integer,intent(inout)::jstack(1:LTE+2) integer::i,j,k,nnpl,nbs,ia,ib,ic,ix real(4)::xa,ya,xb,yb,see nte = 0 nnpl = npl do i = 1,LTE do j = 1,3 ien(i, j) = 0 jee(i, j) = 0 end do end do do if(3 <= nnpl)then nbs = 1 do if(nbs <= nnpl - 1)then ia = nsr(nbs) ib = nsr(nbs + 1) ic = nsr(mod((nbs + 1),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(1.0E-10 < see)then nte = nte + 1 ien(nte, 1) = ia ien(nte, 2) = ib ien(nte, 3) = ic nnpl = nnpl - 1 do i = nbs + 1,nnpl nsr(i) = nsr(i + 1) end do end if nbs = nbs + 1 else exit end if end do else exit end if end do ix = 0 do i = 1,2 * LTE + 1 ihen(i, 1) = 0 ihen(i, 2) = 0 jhen(i) = 0 iad(i) = 0 end do do i = 1,nte loopj: do j = 1,3 ia = ien(i, j) ib = ien(i, mod(j,3) + 1) do k = 1,ix if((ihen(k, 1) == ib).and.(ihen(k, 2) == ia))then jee(i, j) = jhen(k) jee(jhen(k), iad(k)) = i cycle loopj end if end do ix = ix + 1 jhen(ix) = i iad(ix) = j ihen(ix, 1) = ia ihen(ix, 2) = ib end do loopj end do call LAWSON(nte, ien, jee, npl, px, py, jstack) end subroutine SUBDIV !----------------------------------------------------------------------------------------- subroutine LAWSON(nte, ien, jee, npl, px, py, jstack) !----------------------------------------------------------------------------------------- integer,intent(in)::nte integer,intent(inout)::ien(1:LTE,1:3) integer,intent(inout)::jee(1:LTE,1:3) integer,intent(in)::npl real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::jstack(1:LTE+2) integer::i,j,itop,maxstk,ncount,ielm,il,ir,jl1,jl2,jl3,jr1,jr2,jr3 integer::iv1,iv2,iv3,iv4,ia,ib,iswap,iedge real(4)::xx,yy itop = 0 maxstk = npl ncount = 0 do i = 1,nte ielm = i itop = itop + 1 jstack(itop) = IPUSH(ielm, maxstk, itop) end do do if(0 < itop)then ncount = ncount + 1 if(LTE < ncount)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB LAWSON: stop' stop end if il = jstack(itop) itop = itop - 1 do j = 1,3 jl1 = j jl2 = mod(jl1,3) + 1 jl3 = mod(jl2,3) + 1 ir = jee(il, jl1) if(ir == 0)cycle 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 = mod(jr1,3) + 1 jr3 = mod(jr2,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) exit end if end do else exit end if end do end subroutine LAWSON !----------------------------------------------------------------------------------------- subroutine FINE(NODE, px, py, jnb, nei, ifix, NELM, mtj, jac,idm, istack, angl, nedg, kstack, map) !----------------------------------------------------------------------------------------- integer,intent(inout)::NODE real(4),intent(inout)::px(1:KTJ+3) real(4),intent(inout)::py(1:KTJ+3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::ifix(1:KTJ) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer,intent(inout)::istack(1:KTJ) real(4),intent(inout)::angl(1:2*KTJ+1) integer,intent(inout)::nedg(1:2*KTJ+1) integer,intent(inout)::kstack(1:KTJ) integer,intent(inout)::map(1:2*KTJ+1) integer::i,j,k,l,ii,kk,mm real(4)::xx,gx,gy,ar,s,xc,yc,cgrax,cgray,xkp,ykp,drate,almt integer::ntotal,neib,ielm,jelm,kelm,lelm,melm,j1,j2,j3,nn,ls,le,lm,ied,jed,iv1,iv2,jedg ntotal = INPNODT if(ntotal <= NODE)ntotal = NODE almt = 1.0 / real(ntotal) 20 continue almt = 0.5 * almt neib = 0 ! do i = 1,KTE ! angl(i) = 0.0 ! nedg(i) = 0 ! kstack(i) = 0 ! end do do i = 1,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 end do do i = 1,neib - 1 k = kstack(i) xx = angl(k) ii = i + 1 l = i do j = ii,neib kk = kstack(j) if(angl(kk) <= xx)cycle xx = angl(kk) l = j end do mm = kstack(i) kstack(i) = kstack(l) kstack(l) = mm end do 70 continue if(neib == 0)goto 20 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, mod(ied,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 do j = 1,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 end do cgrax = gx / ar cgray = gy / ar xkp = px(NODE) ykp = py(NODE) px(NODE) = cgrax py(NODE) = cgray do j = 1,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 end if end do end if do i = 1,NELM map(i) = 1 end do do i = 1,jnb(NODE) lelm = nei(NODE, i) map(lelm) = 0 angl(lelm) = 0.0 nedg(lelm) = 0 end do nn = 0 do i = 1,neib if(map(kstack(i)) == 1)then nn = nn + 1 kstack(nn) = kstack(i) end if end do neib = nn do i = 1,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 do j = 1,neib kstack(neib - j + 2) = kstack(neib - j + 1) end do neib = neib + 1 kstack(1) = melm else if(drate < angl(kstack(neib)))then neib = neib + 1 kstack(neib) = melm else ls = 1 le = neib do if(le - ls /= 1)then lm = int((ls + le) / 2) if(angl(kstack(lm)) < drate)then le = lm else ls = lm end if else exit end if end do do k = neib,le,-1 kstack(k + 1) = kstack(k) end do neib = neib + 1 kstack(le) = melm end if end if end do goto 70 end if end subroutine FINE !----------------------------------------------------------------------------------------- subroutine REMESH(ielm, ied, jelm, jed, NODE, px, py, jnb, nei, NELM, mtj, jac, idm, istack) !----------------------------------------------------------------------------------------- integer,intent(in)::ielm integer,intent(in)::ied integer,intent(in)::jelm integer,intent(in)::jed integer,intent(in)::NODE real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer,intent(inout)::istack(1:KTJ) integer::itop,maxstk,ip,ibj,jbj,ia,ib,ic,id,iv1,iv2,iv3,iv4,iedge,il,ir,iera,ierb,ierl,iswap real(4)::xp,yp itop = 0 maxstk = NODE ip = NODE xp = px(ip) yp = py(ip) ibj = mod(ied,3) + 1 jbj = mod(jed,3) + 1 ia = jac(ielm, ibj) ib = jac(ielm, mod(ibj,3) + 1) ic = jac(jelm, jbj) id = jac(jelm, mod(jbj,3) + 1) iv1 = mtj(ielm, ibj) iv2 = mtj(ielm, mod(ibj,3) + 1) iv3 = mtj(jelm, jbj) iv4 = mtj(jelm, mod(jbj,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) 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) end if end if if(ic /= 0)then if(idm(ic) == idm(jelm))then itop = itop + 1 istack(itop) = IPUSH(jelm, maxstk, itop) 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) end if end if do if(0 < itop)then il = istack(itop) itop = itop - 1 ir = jac(il, 2) call EDGE(ir, il, jac, ierl) iera = mod(ierl,3) + 1 ierb = mod(iera,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) end if end if if(ib /= 0)then if(idm(ib) == idm(ir))then itop = itop + 1 istack(itop) = IPUSH(ir, maxstk, itop) end if end if if(ic /= 0)then call EDGE(ic, il, jac, iedge) jac(ic, iedge) = ir end if end if else exit end if end do end subroutine REMESH !----------------------------------------------------------------------------------------- subroutine LAPLAS(NODE, ifix, mtj, px, py, jnb, nei) !----------------------------------------------------------------------------------------- integer,intent(in)::NODE integer,intent(in)::ifix(1:KTJ) integer,intent(in)::mtj(1:2*KTJ+1,1:3) real(4),intent(inout)::px(1:KTJ+3) real(4),intent(inout)::py(1:KTJ+3) integer,intent(in)::jnb(1:KTJ) integer,intent(in)::nei(1:KTJ,1:KCM) integer::itera = 5 integer::it,i,j,ielm,j1,j2,j3 real(4)::gx,gy,ar,s,xc,yc,cgrax,cgray do it = 1,itera do i = 1,NODE if(ifix(i) == 0)then gx = 0.0 gy = 0.0 ar = 0.0 do j = 1,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 end do cgrax = gx / ar cgray = gy / ar px(i) = cgrax py(i) = cgray end if end do end do end subroutine LAPLAS !----------------------------------------------------------------------------------------- subroutine REMOVE(NEX, mindex, NELM, mtj, jac, idm) !----------------------------------------------------------------------------------------- integer,intent(in)::NEX integer,intent(inout)::mindex(1:KBD+1) integer,intent(inout)::NELM integer,intent(inout)::mtj(1:2*KTJ+1,1:3) integer,intent(inout)::jac(1:2*KTJ+1,1:3) integer,intent(inout)::idm(1:2*KTJ+1) integer::i,j,k,l,iz,inelm,ielm,jelm,kelm,lelm,ikp,kedg,ledg integer::mkp(1:3),jkp(1:3) ielm = 0 mindex(1) = 1 do i = 1,NEX iz = i inelm = 0 do j = mindex(i),NELM if(idm(j) == iz)then ielm = ielm + 1 jelm = j inelm = inelm + 1 if(ielm /= jelm)then do k = 1,3 mkp(k) = mtj(ielm, k) jkp(k) = jac(ielm, k) end do ikp = idm(ielm) do k = 1,3 kelm = jac(jelm, k) if(kelm /= 0)then call EDGE(kelm, jelm, jac, kedg) jac(kelm, kedg) = ielm + KTE + 1 end if end do do l = 1,3 lelm = jkp(l) if(lelm /= 0)then call EDGE(lelm, ielm, jac, ledg) jac(lelm, ledg) = jelm + KTE + 1 end if end do do k = 1,3 jac(ielm, k) = mod(jac(ielm, k),(KTE + 1)) jac(jelm, k) = mod(jac(jelm, k),(KTE + 1)) kelm = jac(ielm, k) lelm = jac(jelm, k) do l = 1,3 jac(kelm, l) = mod(jac(kelm, l),(KTE + 1)) jac(lelm, l) = mod(jac(lelm, l),(KTE + 1)) end do end do do k = 1,3 jkp(k) = jac(ielm, k) end do do k = 1,3 mtj(ielm, k) = mtj(jelm, k) jac(ielm, k) = jac(jelm, k) mtj(jelm, k) = mkp(k) jac(jelm, k) = jkp(k) end do idm(ielm) = idm(jelm) idm(jelm) = ikp end if end if end do mindex(i + 1) = mindex(i) + inelm end do do i = 1,ielm do j = 1,3 if(ielm < jac(i, j))then jac(i, j) = 0 end if end do end do do i = ielm + 1,NELM do j = 1,3 mtj(i, j) = 0 jac(i, j) = 0 end do end do NELM = ielm end subroutine REMOVE !----------------------------------------------------------------------------------------- subroutine CHECK(NELM, mtj, jac) !----------------------------------------------------------------------------------------- integer,intent(in)::NELM integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer,intent(in)::jac(1:2*KTJ+1,1:3) integer::i,j,k,ielm,jelm,kelm,ia,ib,ja,jb do i = 1,NELM loopj: do j = 1,3 ielm = i ia = mtj(i, j) ib = mtj(i, mod(j,3) + 1) jelm = jac(i, j) if(jelm == 0)cycle do k = 1,3 kelm = jac(jelm, k) if(kelm == ielm)then ja = mtj(jelm, k) jb = mtj(jelm, mod(k,3) + 1) if((ia == jb).and.(ib == ja))cycle loopj !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB CHECK: stop' stop end if end do !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB CHECK: stop' stop end do loopj end do end subroutine CHECK !----------------------------------------------------------------------------------------- subroutine DATA(NODE, NELM, mtj, jac, idm, px, py, fnameW) !----------------------------------------------------------------------------------------- integer,intent(in)::NODE integer,intent(in)::NELM integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer,intent(in)::jac(1:2*KTJ+1,1:3) integer,intent(in)::idm(1:2*KTJ+1) real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) character(len=50),intent(in)::fnameW integer::i,j real(4)::aa open(12, file=fnameW, status='replace') write(12,'(2(i5))') NODE,NELM do i = 1,NELM write(12,'(4(i5))') (mtj(i,j),j=1,3),idm(i) end do do i = 1,NODE write(12,'(2(e15.7),i3)') px(i),py(i) end do ! Find irregular elements do i=1,NELM aa=AREA(i, mtj, px, py) if(aa<=0.0)write(12,'("Irregular element: No, area=",i5,", ",e15.7)'),i,aa end do close(12) end subroutine DATA !----------------------------------------------------------------------------------------- subroutine INCR(n, l, jnb, nei) !----------------------------------------------------------------------------------------- integer,intent(in)::n integer,intent(in)::l integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) if(n <= KTJ)then jnb(n) = jnb(n) + 1 if(KCM < jnb(n))then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB INCR: stop' stop end if nei(n, jnb(n)) = l end if end subroutine INCR !----------------------------------------------------------------------------------------- subroutine DECR(n, l, jnb, nei) !----------------------------------------------------------------------------------------- integer,intent(in)::n integer,intent(in)::l integer,intent(inout)::jnb(1:KTJ) integer,intent(inout)::nei(1:KTJ,1:KCM) integer::i,inb if(n <= KTJ)then inb = NEIBOR(n, l, jnb, nei) do i = inb,jnb(n) - 1 nei(n, i) = nei(n, i + 1) end do nei(n, jnb(n)) = 0 jnb(n) = jnb(n) - 1 if(jnb(n) == 0)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN SUB DECR: stop' stop end if end if end subroutine DECR !----------------------------------------------------------------------------------------- subroutine DISTOR(n, mtj, px, py, drate, jedg) !----------------------------------------------------------------------------------------- integer,intent(in)::n integer,intent(in)::mtj(1:2*KTJ+1,1:3) real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) real(4),intent(inout)::drate integer,intent(inout)::jedg real(4)::dst(1:3) real(4)::xa,ya,xb,yb,xc,yc,ang1,ang2,ang3 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 = abs(dst(1)) + abs(dst(2)) + abs(dst(3)) jedg = 1 if(dst(jedg) < dst(2))jedg = 2 if(dst(jedg) < dst(3))jedg = 3 end subroutine DISTOR !----------------------------------------------------------------------------------------- subroutine LOCATE(xp, yp, px, py, mtj, jac, NELM, itri) !----------------------------------------------------------------------------------------- real(4),intent(in)::xp real(4),intent(in)::yp real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer,intent(in)::jac(1:2*KTJ+1,1:3) integer,intent(in)::NELM integer,intent(inout)::itri integer::i,ia,ib itri = NELM 10 continue do i = 1,3 ia = mtj(itri, i) ib = mtj(itri, mod(i,3) + 1) if((px(ia) - xp) * (py(ib) - yp) < (py(ia) - yp) * (px(ib) - xp))then itri = jac(itri, i) goto 10 end if end do end subroutine LOCATE !----------------------------------------------------------------------------------------- subroutine SWAP(x1, y1, x2, y2, x3, y3, xp, yp, iswap) !----------------------------------------------------------------------------------------- real(4),intent(in)::x1 real(4),intent(in)::y1 real(4),intent(in)::x2 real(4),intent(in)::y2 real(4),intent(in)::x3 real(4),intent(in)::y3 real(4),intent(in)::xp real(4),intent(in)::yp integer,intent(inout)::iswap real(4)::x13,y13,x23,y23,x1p,y1p,x2p,y2p,cosa,cosb,sina,sinb 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 else if((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.0)then iswap = 1 else iswap = 0 end if end if end subroutine SWAP !----------------------------------------------------------------------------------------- subroutine EDGE(l, k, work, iedge) !----------------------------------------------------------------------------------------- integer,intent(in)::l integer,intent(in)::k integer,intent(in)::work(:,:) integer,intent(inout)::iedge integer::i iedge = -1 do i = 1,3 if(work(l, i) == k)then iedge = i exit end if end do !Show ERROR mesage and stop the execution if(iedge == -1)then write(6,*) 'ERROR IN SUB EDGE: stop' stop end if end subroutine EDGE !----------------------------------------------------------------------------------------- integer function IPUSH(item, maxstk, itop) !----------------------------------------------------------------------------------------- integer,intent(in)::item integer,intent(in)::maxstk integer,intent(in)::itop if(maxstk < itop)then !Show ERROR mesage and stop the execution write(6,*) 'ERROR IN FUNCTION IPUSH: stop' stop else IPUSH = item end if end function IPUSH !----------------------------------------------------------------------------------------- integer function NEIBOR(n, l, jnb, nei) !----------------------------------------------------------------------------------------- integer,intent(in)::n integer,intent(in)::l integer,intent(in)::jnb(1:KTJ) integer,intent(in)::nei(1:KTJ,1:KCM) integer::i NEIBOR = -1 do i = 1,jnb(n) if(nei(n, i) == l)then NEIBOR = i exit end if end do !Show ERROR mesage and stop the execution if(NEIBOR == -1)then write(6,*) 'ERROR IN FUNCTION NEIBOR: stop' stop end if end function NEIBOR !----------------------------------------------------------------------------------------- integer function IVERT(l, k, mtj) !----------------------------------------------------------------------------------------- integer,intent(in)::l integer,intent(in)::k integer,intent(in)::mtj(1:2*KTJ+1,1:3) integer::i IVERT = -1 do i = 1,3 if(mtj(l, i) == k)then IVERT = i exit end if end do !Show ERROR mesage and stop the execution if(IVERT == -1)then write(6,*) 'ERROR IN FUNCTION IVERT: stop' stop end if end function IVERT !----------------------------------------------------------------------------------------- real(4) function AREA(ielm, mtj, px, py) !----------------------------------------------------------------------------------------- integer,intent(in)::ielm integer,intent(in)::mtj(1:2*KTJ+1,1:3) real(4),intent(in)::px(1:KTJ+3) real(4),intent(in)::py(1:KTJ+3) integer::j1,j2,j3 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 AREA !----------------------------------------------------------------------------------------- real(4) function THETA(x0, y0, x1, y1, x2, y2) !----------------------------------------------------------------------------------------- real(4),intent(in)::x0 real(4),intent(in)::y0 real(4),intent(in)::x1 real(4),intent(in)::y1 real(4),intent(in)::x2 real(4),intent(in)::y2 real(4)::err,xa,ya,xb,yb,prdin,prdex err = 1.0E-10 xa = x1 - x0 ya = y1 - y0 xb = x2 - x0 yb = y2 - y0 prdin = xa * xb + ya * yb prdex = xa * yb - xb * ya if(abs(prdin) < err)then THETA = pi / 2.0 else THETA = 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 THETA !----------------------------------------------------------------------------------------- end program f90_FEM_MESH3 !-----------------------------------------------------------------------------------------