module defparam_C
    implicit none
    integer,parameter::KTE = 20000
    integer,parameter::KTJ = 20000
    integer,parameter::KED = 20000
    integer,parameter::KCM = 100
end module defparam_C

!-----------------------------------------------------------------------------------------
program f90_fem_mesh4c
!-----------------------------------------------------------------------------------------
    use defparam_C
    implicit none
    integer::NODE           =0
    integer::NELM           =0
    integer::mtj(1:KTE,1:4) =0
    real(8)::px(1:KTJ)      =0.0D0
    real(8)::py(1:KTJ)      =0.0D0
    integer::ifix(1:KTJ)    =0
    character(len=50)::fnameR
    character(len=50)::fnameW

    call getarg(1,fnameR)
    call getarg(2,fnameW)
    call SQINPUT(NODE, NELM, mtj, px, py, ifix, fnameR)
    call SQLAPLAS(NODE, NELM, mtj, px, py, ifix)
    call SQDATA(NODE, NELM, mtj, px, py, fnameW)

stop

contains
!-----------------------------------------------------------------------------------------
    subroutine SQINPUT(NODE, NELM, mtj, px, py, ifix, fnameR)
!-----------------------------------------------------------------------------------------
        integer,intent(inout)::NODE
        integer,intent(inout)::NELM
        integer,intent(inout)::mtj(1:KTE,1:4)
        real(8),intent(inout)::px(1:KTJ)
        real(8),intent(inout)::py(1:KTJ)
        integer,intent(inout)::ifix(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,4)
        end do
        do i = 1,NODE
            read(11,*) px(i),py(i),ifix(i)
        end do
        close(11)

    end subroutine SQINPUT
!-----------------------------------------------------------------------------------------
    subroutine SQLAPLAS(NODE, NELM, mtj, px, py, ifix)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NODE
        integer,intent(in)::NELM
        integer,intent(inout)::mtj(1:KTE,1:4)
        real(8),intent(inout)::px(1:KTJ)
        real(8),intent(inout)::py(1:KTJ)
        integer,intent(inout)::ifix(1:KTJ)

        integer::i,j,k,l,ja,jb,ip,ix,ka,kb,la,lb,itera,it,j1,j2,j3,j4,ielm
        integer::iflag,iarea
        real(8)::gx,gy,ar,s,xc,yc,cgrax,cgray
        integer::jac(1:KED,1:4)  =0
        integer::ihen(1:KED,1:2) =0
        integer::jhen(1:KED)     =0
        integer::khen(1:KED)     =0
        integer::jnb(1:KED)      =0
        integer::nei(1:KED,1:KCM)=0

        ix=0
        do j = 1,NELM
            do k = 1,4
                ka = mtj(j, k)
                kb = mtj(j, mod(k,4) + 1)
                iflag=0
                ix = ix + 1
                do l = 1,ix
                    la = ihen(l, 1)
                    lb = ihen(l, 2)
                    if((ka == lb).and.(kb == la))then
                        iflag=1
                        exit
                    end if
                end do
                if(iflag==0)then
                    ihen(ix, 1) = ka
                    ihen(ix, 2) = kb
                    jhen(ix) = j
                    khen(ix) = k
                else
                    jac(j, k) = jhen(l)
                    jac(jhen(l), khen(l)) = j
                end if
            end do
        end do

        do i = 1,NELM
            do j = 1,4
                if(jac(i, j) == 0)then
                    ja = mtj(i, j)
                    jb = mtj(i, mod(j,4) + 1)
                    ifix(ja) = 1
                    ifix(jb) = 1
                end if
            end do
        end do

        do i = 1,NELM
            do j = 1,4
                ip = mtj(i, j)
                jnb(ip) = jnb(ip) + 1
                nei(ip, jnb(ip)) = i
            end do
        end do

        itera = 5
        do it = 1,itera
            do i = 1,NODE
                if(ifix(i) == 0)then
                    gx = 0.0D0
                    gy = 0.0D0
                    ar = 0.0D0
                    do j = 1,jnb(i)
                        ielm = nei(i, j)
                        j1 = mtj(ielm, 1)
                        j2 = mtj(ielm, 2)
                        j3 = mtj(ielm, 3)
                        j4 = mtj(ielm, 4)
                        s = SQAREA(ielm, mtj, px, py, iarea)
                        xc = (px(j1) + px(j2) + px(j3) + px(j4)) / 4.0
                        yc = (py(j1) + py(j2) + py(j3) + py(j4)) / 4.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 SQLAPLAS
!-----------------------------------------------------------------------------------------
    subroutine SQDATA(NODE, NELM, mtj, px, py, fnameW)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::NODE
        integer,intent(in)::NELM
        integer,intent(in)::mtj(1:KTE,1:4)
        real(8),intent(in)::px(1:KTJ)
        real(8),intent(in)::py(1:KTJ)
        character(len=50),intent(in)::fnameW

        integer::i,j,iarea
        real(8)::area

        open(12, file=fnameW, status='replace')
        write(12,'(2(i5))') NODE,NELM
        do i = 1,NELM
            write(12,'(7(i5))') (mtj(i,j),j=1,4)
        end do
        do i = 1,NODE
            write(12,'(2(e15.7))') px(i),py(i)
        end do
        ! Find irregular elements
        do i=1,NELM
            area=SQAREA(i, mtj, px, py, iarea)
            if(iarea==1)then
                write(12,'("Irregular element: No,area,iarea=",i5,", ",e15.7,", ",i5)'),i,area,iarea
                do j=1,4
                write(12,'(i5,2(e15.7))'),mtj(i,j),px(mtj(i,j)),py(mtj(i,j))
                end do
            end if
        end do
        close(12)

    end subroutine SQDATA
!-----------------------------------------------------------------------------------------
    real(8) function SQAREA(ielm, mtj, px, py, iarea)
!-----------------------------------------------------------------------------------------
        integer,intent(in)::ielm
        integer,intent(in)::mtj(1:KTE,1:4)
        real(8),intent(in)::px(1:KTJ)
        real(8),intent(in)::py(1:KTJ)
        integer,intent(out)::iarea

        integer::j1,j2,j3,j4
        real(8)::area1,area2,area3,area4

        iarea=1
        j1 = mtj(ielm, 1)
        j2 = mtj(ielm, 2)
        j3 = mtj(ielm, 3)
        j4 = mtj(ielm, 4)
        area1 = 0.5D0 * &
        & (px(j1) * py(j2) + px(j2) * py(j3) + px(j3) * py(j1) - px(j1) * py(j3) - px(j2) * py(j1) - px(j3) * py(j2))
        area2 = 0.5D0 * &
        & (px(j1) * py(j3) + px(j3) * py(j4) + px(j4) * py(j1) - px(j1) * py(j4) - px(j3) * py(j1) - px(j4) * py(j3))
        j1 = mtj(ielm, 2)
        j2 = mtj(ielm, 3)
        j3 = mtj(ielm, 4)
        j4 = mtj(ielm, 1)
        area3 = 0.5D0 * &
        & (px(j1) * py(j2) + px(j2) * py(j3) + px(j3) * py(j1) - px(j1) * py(j3) - px(j2) * py(j1) - px(j3) * py(j2))
        area4 = 0.5D0 * &
        & (px(j1) * py(j3) + px(j3) * py(j4) + px(j4) * py(j1) - px(j1) * py(j4) - px(j3) * py(j1) - px(j4) * py(j3))
        if((0.0D0<area1).and.(0.0D0<area2).and.(0.0D0<area3).and.(0.0D0<area4))iarea=0
        SQAREA = area1 + area2
    end function SQAREA
!-----------------------------------------------------------------------------------------
end program f90_fem_mesh4c
!-----------------------------------------------------------------------------------------
