program f90_fem_domain
    implicit none
    integer,parameter::ntmax=10
    integer,parameter::mdmax=100
    integer::NODT,NELT,MATEL,nod
    integer,allocatable::kakom(:,:)
    integer,allocatable::matno(:)
    real(8),allocatable::x(:),y(:)
    real(8),allocatable::xg(:),yg(:)
    integer::nt(1:ntmax)
    real(8):: x_dom(1:ntmax,1:mdmax),y_dom(1:ntmax,1:mdmax)
    integer::i,j,ne
    real(8)::sumx,sumy
    character(len=50)::fnameR1
    character(len=50)::fnameR2
    character(len=50)::fnameW
    character(len=50)::dummy

    call getarg(1,dummy)
    call getarg(2,fnameR1)
    call getarg(3,fnameR2)
    call getarg(4,fnameW)
    read(dummy,*) nod

    ! Input of Mesh data
    open(11,file=fnameR1,status='old')
    read(11,*) NODT,NELT
    allocate(kakom(1:NELT,1:nod),matno(1:NELT))
    allocate(x(1:NODT),y(1:NODT))
    allocate(xg(1:NELT),yg(1:NELT))
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,nod)
    end do
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    close(11)

    do ne=1,NELT
        sumx=0.0D0
        sumy=0.0D0
        do i=1,nod
            sumx=sumx+x(kakom(ne,i))
            sumy=sumy+y(kakom(ne,i))
        end do
        xg(ne)=sumx/dble(nod)
        yg(ne)=sumy/dble(nod)
    end do

    ! Input of Domain data
    open(11,file=fnameR2,status='old')
    read(11,*) MATEL
    do i=1,MATEL
        read(11,*) nt(i)
        do j=1,nt(i)
            read(11,*) x_dom(i,j),y_dom(i,j)
        end do
    end do
    close(11)

    ! Definition of Domain (material No.)
    call DEFDOM(NELT,MATEL,nt,x_dom,y_dom,xg,yg,matno,ntmax,mdmax)

    open(12, file=fnameW, status='replace')
    write(12,'(4(i6))') nod,NODT,NELT,MATEL
    if(nod==3)then
        do ne=1,NELT
            write(12,'(4(i6))') (kakom(ne,j),j=1,nod),matno(ne)
        end do
    end if
    if(nod==4)then
        do ne=1,NELT
            write(12,'(5(i6))') (kakom(ne,j),j=1,nod),matno(ne)
        end do
    end if
    do i=1,NODT
        write(12,'(2(e15.7))') x(i),y(i)
    end do
    close(12)

    do i=1,MATEL
        dummy(1:50)=''
        write (dummy,'("_mat_",i2.2,".txt")') i
        open(12,file=dummy,status='replace')
        do ne=1,NELT
            if(matno(ne)==i)then
                write(12,'(a)') '>'
                do j=1,nod
                    write(12,'(2(e15.7))') x(kakom(ne,j)),y(kakom(ne,j))
                end do
            end if
        end do
        close(12)
    end do

stop

contains

    subroutine DEFDOM(NELT,MATEL,nt,x_dom,y_dom,xg,yg,matno,ntmax,mdmax)
        integer,intent(in)::NELT,MATEL,ntmax,mdmax
        integer,intent(in)::nt(1:ntmax)
        real(8),intent(in):: x_dom(1:ntmax,1:mdmax),y_dom(1:ntmax,1:mdmax)
        real(8),intent(in)::xg(1:NELT),yg(1:NELT)
        integer,intent(out)::matno(1:NELT)

        real(8)::theta,dx1,dy1,dx2,dy2,r1,r2,cs,sn
        integer::i,j,ne,j1,j2
        real(8),parameter::pi=3.14159265358979323846D0
        real(8),parameter::eps=1D-10

        do ne=1,NELT
            matno(ne)=0
            do i=1,MATEL
                theta=0.0D0
                do j=1,nt(i)
                    j1=j
                    j2=j+1
                    if(j1==nt(i))j2=1
                    dx1=x_dom(i,j1)-xg(ne)
                    dy1=y_dom(i,j1)-yg(ne)
                    dx2=x_dom(i,j2)-xg(ne)
                    dy2=y_dom(i,j2)-yg(ne)
                    r1=sqrt(dx1*dx1+dy1*dy1)
                    r2=sqrt(dx2*dx2+dy2*dy2)
                    if(eps<r1.and.eps<r2)then
                        cs=(dx1*dx2+dy1*dy2)/r1/r2 !inner product
                        sn=(dx1*dy2-dx2*dy1)/r1/r2 !outer product
                        theta=theta+atan2(sn,cs)/pi*180.0D0
                    else
                        theta=theta
                    end if
                end do
                if(eps<abs(theta))matno(ne)=i
            end do
        end do
    end subroutine DEFDOM

end program f90_fem_domain
