module defparam_B
    implicit none
    integer,parameter::KTE = 10000
    integer,parameter::KTJ = 10000
    integer,parameter::KED = 10000
    integer,parameter::KCM = 100
end module defparam_B

!-----------------------------------------------------------------------------------------
program f90_FEM_MESH4B
!-----------------------------------------------------------------------------------------
    use defparam_B
    implicit none
    integer::NODE            =0
    integer::NELM            =0
    integer::mtj(1:KTE,1:9)  =0
    real(8)::px(1:KTJ)       =0.0D0
    real(8)::py(1:KTJ)       =0.0D0
    integer::ifix(1:KTJ)     =0
    integer::jac(1:KTE,1:4)  =0
    integer::idm(1:KTE)      =0
    integer::mmtj(1:KTE,1:4) =0
    integer::kv(1:KTE)       =0
    integer::id(1:KTE)       =0
    integer::map(1:KTE)      =0
    character(len=50)::fnameR
    character(len=50)::fnameW

    call getarg(1,fnameR)
    call getarg(2,fnameW)
    call QUINPUT(NODE, NELM, mtj, jac, idm, ifix, px, py, fnameR)
    call QUEXELM(NELM, mtj, jac, kv, idm, map, px, py)
    call QUISOGEN(NODE, NELM, mtj, jac, px, py, id, mmtj, ifix, idm)
    call QUDATA(NODE, NELM, mmtj, px, py, ifix, fnameW)

stop

contains
!-----------------------------------------------------------------------------------------
    subroutine QUINPUT(NODE, NELM, mtj, jac, idm, ifix, px, py, fnameR)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NODE
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:9)
        integer,intent(inout)::jac(1:KTE,1:4)
        integer,intent(inout)::idm(1:KTE)
        integer,intent(inout)::ifix(1:KTJ)
        real(8),intent(inout)::px(1:KTJ)
        real(8),intent(inout)::py(1:KTJ)
        character(len=50),intent(in)::fnameR

        integer::i,j

        open(11,file=fnameR,status='old')
        read(11,*) NODE,NELM
        do i = 1,NELM
            read(11,*) (mtj(i,j),j=1,3),(jac(i,j),j=1,3),idm(i)
        end do
        do i = 1,NODE
            read(11,*) px(i),py(i),ifix(i)
        end do
        close(11)

    end subroutine QUINPUT
!-----------------------------------------------------------------------------------------
    subroutine QUEXELM(NELM, mtj, jac, kv, idm, map, px, py)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:9)
        integer,intent(inout)::jac(1:KTE,1:4)
        integer,intent(inout)::kv(1:KTE)
        integer,intent(inout)::idm(1:KTE)
        integer,intent(inout)::map(1:KTE)
        real(8),intent(inout)::px(1:KTJ)
        real(8),intent(inout)::py(1:KTJ)

        integer::i,j,k,iv,ip,jn,ie1,ie2,ie3,ia,ib,ic,jacia,jacib,inn1,inn2,ien1,ien2
        real(8)::pxa,pya,pxb,pyb,pxc,pyc,px12,py12,px23,py23,pl12,pl23
        real(8)::pl(1:3)

        iv = 0
        do i = 1,NELM
            if(mtj(i, 1) == 0)cycle
            if(mtj(i, 4) /=0)cycle
            pxa = px(mtj(i, 1)) - px(mtj(i, 2))
            pya = py(mtj(i, 1)) - py(mtj(i, 2))
            pxb = px(mtj(i, 2)) - px(mtj(i, 3))
            pyb = py(mtj(i, 2)) - py(mtj(i, 3))
            pxc = px(mtj(i, 3)) - px(mtj(i, 1))
            pyc = py(mtj(i, 3)) - py(mtj(i, 1))
            pl(1) = pxa * pxa + pya * pya
            pl(2) = pxb * pxb + pyb * pyb
            pl(3) = pxc * pxc + pyc * pyc
            do j = 1,3
                if((pl(mod(j,3) + 1) < pl(j)).and.(pl(mod((j + 1),3) + 1) < pl(j)))then
                    ip = j
                    exit
                end if
            end do
            jn = jac(i, ip)
            if(idm(i) /=idm(jn))cycle
            if(mtj(jn, 4) /=0)cycle
            if(jn == 0)cycle
            call QUEDGE(jn, i, jac, ie1)
            ie2 = mod(ie1,3) + 1
            ie3 = mod(ie2,3) + 1
            px12 = px(mtj(jn, ie2)) - px(mtj(jn, ie3))
            py12 = py(mtj(jn, ie2)) - py(mtj(jn, ie3))
            px23 = px(mtj(jn, ie3)) - px(mtj(jn, ie1))
            py23 = py(mtj(jn, ie3)) - py(mtj(jn, ie1))
            pl12 = px12 * px12 + py12 * py12
            pl23 = px23 * px23 + py23 * py23
            if((pl(ip) < pl12).or.(pl(ip) < pl23))cycle

            ia = mtj(i, mod(ip,3) + 1)
            ib = mtj(i, mod((ip + 1),3) + 1)
            ic = mtj(i, ip)
            jacia = jac(i, mod(ip,3) + 1)
            jacib = jac(i, mod((ip + 1),3) + 1)
            mtj(i, 1) = ia
            mtj(i, 2) = ib
            mtj(i, 3) = ic
            mtj(i, 4) = mtj(jn, ie3)
            jac(i, 1) = jacia
            jac(i, 2) = jacib
            inn1 = jac(jn, ie2)
            inn2 = jac(jn, ie3)

            if(inn1 /=0)then
                 call QUEDGE(inn1, jn, jac, ien1)
                jac(inn1, ien1) = i
            end if

            if(inn2 /= 0)then
                call QUEDGE(inn2, jn, jac, ien2)
                jac(inn2, ien2) = i
            end if
            jac(i, 3) = inn1
            jac(i, 4) = inn2

            do k = 1,3
                mtj(jn, k) = 0
                jac(jn, k) = 0
                idm(jn) = 0
            end do
            iv = iv + 1
            kv(iv) = jn
        end do

         call QUDELETE(NELM, mtj, jac, idm, iv, kv, map)

    end subroutine QUEXELM
!-----------------------------------------------------------------------------------------
    subroutine QUDELETE(NELM, mtj, jac, idm, iv, kv, map)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:9)
        integer,intent(inout)::jac(1:KTE,1:4)
        integer,intent(inout)::idm(1:KTE)
        integer,intent(inout)::iv
        integer,intent(inout)::kv(1:KTE)
        integer,intent(inout)::map(1:KTE)

        integer::i,m,n,ia

        m = 0
        n = 0
        do i = 1,NELM
            map(i) = 1
        end do
        do i = 1,iv
            map(kv(i)) = 0
        end do
        do i = 1,NELM
            if(map(i) /=0)then
                m = m + 1
                map(i) = m
            end if
        end do

        do i = 1,NELM
            if(map(i) /=0)then
                n = n + 1
                do ia = 1,4
                    mtj(n, ia) = mtj(i, ia)
                    idm(n) = idm(i)
                    if(jac(i, ia) == 0)then
                        jac(n, ia) = 0
                    else
                        jac(n, ia) = map(jac(i, ia))
                    end if
                end do
            end if
        end do

        do i = n + 1,NELM
            do ia = 1,4
                mtj(i, ia) = 0
                jac(i, ia) = 0
                idm(i) = 0
            end do
        end do

        NELM = NELM - iv

    end subroutine QUDELETE
!-----------------------------------------------------------------------------------------
    subroutine QUISOGEN(NODE, NELM, mtj, jac, px, py, id, mmtj, ifix, idm)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NODE
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:9)
        integer,intent(inout)::jac(1:KTE,1:4)
        real(8),intent(inout)::px(1:KTJ)
        real(8),intent(inout)::py(1:KTJ)
        integer,intent(inout)::id(1:KTE)
        integer,intent(inout)::mmtj(1:KTE,1:4)
        integer,intent(inout)::ifix(1:KTJ)
        integer,intent(inout)::idm(1:KTE)

        integer::i,j,k,jn,ien
        real(8)::pxx,pyy

        do i = 1,NELM
            if(mtj(i, 4) == 0)then
                id(i) = 3
                mtj(i, 5) = mtj(i, 3)
                mtj(i, 3) = mtj(i, 2)
                mtj(i, 2) = 0
            else
                id(i) = 4
                mtj(i, 7) = mtj(i, 4)
                mtj(i, 5) = mtj(i, 3)
                mtj(i, 4) = 0
                mtj(i, 3) = mtj(i, 2)
                mtj(i, 2) = 0
            end if
        end do

        do i = 1,NELM
            do j = 1,id(i)
                if(mtj(i, 2 * j) /=0)cycle
                NODE = NODE + 1
                mtj(i, 2 * j) = NODE
                px(NODE) = (px(mtj(i, 2 * j - 1)) + px(mtj(i, mod((2 * j),(2 * id(i))) + 1))) / 2.0D0
                py(NODE) = (py(mtj(i, 2 * j - 1)) + py(mtj(i, mod((2 * j),(2 * id(i))) + 1))) / 2.0D0
                if(jac(i, j) == 0)then
                    ifix(NODE) = 1
                    cycle
                end if
                jn = jac(i, j)
                call QUEDGE(jn, i, jac, ien)
                if(mtj(jn, 2 * ien) /=0)cycle
                mtj(jn, 2 * ien) = NODE
                if(idm(i) /=idm(jn))ifix(NODE) = 1
            end do

            NODE = NODE + 1
            mtj(i, 2 * id(i) + 1) = NODE
            pxx = 0.0D0
            pyy = 0.0D0
            do k = 1,id(i)
                pxx = pxx + px(mtj(i, 2 * k))
                pyy = pyy + py(mtj(i, 2 * k))
            end do
            px(NODE) = pxx / real(id(i))
            py(NODE) = pyy / real(id(i))
        end do

         call QUSQUARE(NELM, mtj, id, mmtj)

    end subroutine QUISOGEN
!-----------------------------------------------------------------------------------------
    subroutine QUEDGE(l, k, jac, iedge)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::l
        integer,intent(in)::k
        integer,intent(in)::jac(1:KTE,1:4)
        integer,intent(out)::iedge

        integer::i

        iedge = -1
        do i = 1,4
            if(jac(l, i) == k)then
                iedge = i
                exit
            end if
        end do
        !ERRORメッセージ＆実行中止
        if(iedge == -1)then
            write(6,*) 'ERROR IN SUB QUEDGE: stop'
            stop
        end if
    end subroutine QUEDGE
!-----------------------------------------------------------------------------------------
    subroutine QUSQUARE(NELM, mtj, id, mmtj)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:9)
        integer,intent(inout)::id(1:KTE)
        integer,intent(inout)::mmtj(1:KTE,1:4)

        integer::i,j,nelm1

        nelm1 = 0
        do i = 1,NELM
            if(mtj(i, 9) == 0)then
                do j = 1,id(i)
                    nelm1 = nelm1 + 1
                    mmtj(nelm1, 1) = mtj(i, mod((2 * j + 3),6) + 1)
                    mmtj(nelm1, 2) = mtj(i, 2 * j - 1)
                    mmtj(nelm1, 3) = mtj(i, 2 * j)
                    mmtj(nelm1, 4) = mtj(i, 7)
                end do
            else
                do j = 1,id(i)
                    nelm1 = nelm1 + 1
                    mmtj(nelm1, 1) = mtj(i, mod((2 * j + 5),8) + 1)
                    mmtj(nelm1, 2) = mtj(i, 2 * j - 1)
                    mmtj(nelm1, 3) = mtj(i, 2 * j)
                    mmtj(nelm1, 4) = mtj(i, 9)
                end do
            end if
        end do

        NELM = nelm1

    end subroutine QUSQUARE
!-----------------------------------------------------------------------------------------
    subroutine QUDATA(NODE, NELM, mmtj, px, py, ifix, fnameW)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NODE
        integer,intent(in)::NELM
        integer,intent(in)::mmtj(1:KTE,1:4)
        real(8),intent(in)::px(1:KTJ)
        real(8),intent(in)::py(1:KTJ)
        integer,intent(in)::ifix(1:KTJ)
        character(len=50),intent(in)::fnameW

        integer::i,j

        open(12, file=fnameW, status='replace')
        write(12,'(2(i5))') NODE,NELM
        do i = 1,NELM
            write(12,'(4(i5))') (mmtj(i,j),j=1,4)
        end do
        do i = 1,NODE
            write(12,'(2(e15.7),i3)') px(i),py(i),ifix(i)
        end do
        close(12)

    end subroutine QUDATA
!-----------------------------------------------------------------------------------------
end program f90_FEM_MESH4B
!-----------------------------------------------------------------------------------------