!----------------------------------------------------------------------
! 方形導波管のモード関数を有限要素法で計算するプログラム
!                                                           By 平野拓一
!
! 言語: FORTRAN90
! コンパイル環境: Compaq Visual Fortran
!----------------------------------------------------------------------

!======================================================================
! 定数
!======================================================================
module mod_consts
   real(8),parameter :: pi=3.141592653589793d0       ! 円周率
   real(8),parameter :: velocity_of_light=2.998d8    ! 光速
   real(8),parameter :: z0=120.0d0*pi
   real(8),parameter :: y0=1.0d0/(120.0d0*pi)
   complex(8),parameter :: zj=dcmplx(0.0d0,1.0d0)
end module

!======================================================================
! FEM
!======================================================================
module mod_fem
    !-------- 要素に関する情報のデータ型 --------
    type element
        integer :: node(3)    ! ローカルノード → グローバルノード
    end type element

    ! 配列を固定でとるために最大数を指定する
    integer,parameter :: mx=20                      ! x 方向最大分割数
    integer,parameter :: my=10                      ! y 方向最大分割数
    integer,parameter :: melem=2*mx*my              ! 最大要素数
    integer,parameter :: mnode=(mx+1)*(my+1)        ! 最大ノード数

    type (element) :: elem(melem)
    real(8) :: gn2x(mnode),gn2y(mnode)
    real(8) :: sizex,sizey
    integer :: nx,ny,nelem,nnode,npec,npmc

    real(8) :: syst(mnode,mnode),sysu(mnode,mnode)

    real(8) :: eigval(mnode),eigvec(mnode,mnode)    ! 固有値問題の解
end module

!**********************************************************************
! メインプログラム
!**********************************************************************

program rect
    use mod_fem
    implicit none

    call input                 ! 入力ファイル読みこみ
    call mesh                  ! メッシュデータを作成
    call testmesh              ! メッシュデータの確認ファイル出力
    call sysmatrix             ! 系行列を作成
    call eigsolv               ! 固有値問題を解く
    call outeig                ! 固有値をファイルに出力
    call outpot(2)             ! モード関数（ポテンシャルデータ）をファイルに出力
    call outfield(2)           ! モード関数（フィールドデータ）をファイルに出力

    stop
end program

!----------------------------------------------------------------------
! ファイル入力
!----------------------------------------------------------------------
subroutine input
   use mod_fem
   implicit none

   integer :: i
   character(128) :: text

   write(*,*) '**** INPUT ****'

   open(10,file='input.dat')
   read(10,*) text,nx
   read(10,*) text,ny
   read(10,*) text,sizex
   read(10,*) text,sizey
   close(10)

   return
end subroutine

!----------------------------------------------------------------------
! メッシュデータの作成
!
! [Input]
! nx: x 方向の分割数
! ny: y 方向の分割数
! sizex: x 方向の長さ [m]
! sizey: y 方向の長さ [m]
!
! [Output]
! nelem: 要素数
! nnode: 系全体の節点数
! nx: x 方向の分割数
! ny: y 方向の分割数
! cut: メッシュの切り方 (0 or 1)
! l2gid: 要素 ie の局所節点番号 → 広域節点番号
! x: 広域節点番号の x 座標
! y: 広域節点番号の y 座標
!----------------------------------------------------------------------
subroutine mesh
    use mod_fem
    implicit none

    integer :: i,j,n,k1,k2,k3,k4

    write(*,*) '**** MESH ****'

    nelem=nx*ny*2             ! 三角要素の数
    nnode=(nx+1)*(ny+1)       ! 系全体の節点数

    ! ローカル→グローバルへの変換
    !
    ! k4                              k3
    ! **********************************
    ! *                            *****
    ! *                        ****    *
    ! *                    ****        *
    ! *                ****            *
    ! *            ****                *
    ! *        ****                    *
    ! *    ****                        *
    ! *****                            *
    ! **********************************
    ! k1                              k2
    !
    n=0
    do i=1,nx
        do j=1,ny
            k1=(ny+1)*(i-1)+(j  )    ! 左下
            k2=(ny+1)*(i  )+(j  )    ! 右下
            k3=(ny+1)*(i  )+(j+1)    ! 右上
            k4=(ny+1)*(i-1)+(j+1)    ! 左上

            n=n+1
            elem(n).node(1)=k1
            elem(n).node(2)=k2
            elem(n).node(3)=k3
            n=n+1
            elem(n).node(1)=k3
            elem(n).node(2)=k4
            elem(n).node(3)=k1  
        end do
    end do

    ! 節点の x,y 座標を格納する
    n=0
    do i=1,nx+1
        do j=1,ny+1
            n=n+1
            gn2x(n)=sizex*dfloat(i-1)/dfloat(nx)
            gn2y(n)=sizey*dfloat(j-1)/dfloat(ny)
        end do
    end do

    return
end subroutine

!----------------------------------------------------------------------
! 要素行列の計算
!
! [Input]
! ie: 対象とする要素の要素番号
!
! [Output]
! t: K∇
! u: K
!----------------------------------------------------------------------
subroutine elemmatrix(ie,t,u)
    use mod_fem
    implicit none

    integer :: ie
    real(8) :: t(3,3),u(3,3)

    integer :: i,j,n1,n2,n3
    real(8) :: x(3),y(3),b(3),c(3),area

    n1=elem(ie).node(1)
    n2=elem(ie).node(2)
    n3=elem(ie).node(3)

    x(1)=gn2x(n1)
    x(2)=gn2x(n2)
    x(3)=gn2x(n3)
    y(1)=gn2y(n1)
    y(2)=gn2y(n2)
    y(3)=gn2y(n3)

    b(1)=y(2)-y(3)
    b(2)=y(3)-y(1)
    b(3)=y(1)-y(2)
    c(1)=x(3)-x(2)
    c(2)=x(1)-x(3)
    c(3)=x(2)-x(1)
    area=(x(1)*b(1)+x(2)*b(2)+x(3)*b(3))/2.0d0

    ! エラー処理
    if(area.le.0.0d0) then
       write(*,*) 'Sub(elemmatrix):'
       write(*,*) '  area <= 0'
       write(*,*) '  element number=',ie
       stop
    endif

    do i=1,3
        do j=1,3
            ! T=K∇ 行列の計算
            t(i,j)=(b(i)*b(j)+c(i)*c(j))/(area*4.0d0)

            ! U=K 行列の計算
            if(i.eq.j) then
               u(i,j)=area/6.0d0
            else
               u(i,j)=area/12.0d0
            end if
        end do
    end do

    return
end subroutine

!----------------------------------------------------------------------
! 系行列の作成
!
! [Output]
! mt: K∇
! mu: K
!----------------------------------------------------------------------
subroutine sysmatrix
    use mod_fem
    implicit none

    real(8) :: t(3,3),u(3,3)
    integer :: i,j,ie,n(3)

    write(*,*) '**** SYSMATRIX ****'

    do i=1,nnode
        do j=1,nnode
            syst(i,j)=0.0d0
            sysu(i,j)=0.0d0
        end do
    end do

    do ie=1,nelem
        call elemmatrix(ie,t,u)
        n(1)=elem(ie).node(1)
        n(2)=elem(ie).node(2)
        n(3)=elem(ie).node(3)
        do i=1,3
            do j=1,3
                syst(n(i),n(j))=syst(n(i),n(j))+t(i,j)
                sysu(n(i),n(j))=sysu(n(i),n(j))+u(i,j)
            end do
        end do
    end do

    return
end subroutine

!----------------------------------------------------------------------
! 固有値問題を解く
!----------------------------------------------------------------------
subroutine eigsolv
    use mod_fem
    implicit none

    integer,parameter :: nn1=mnode
    integer :: nne,nnv
    real(8) :: a(nn1,nn1),b(nn1,nn1),v(nn1,nn1),e(nn1)
    real(8) :: w(mnode,7)
    integer :: i,j,ier

    write(*,*) '**** SOLVE ****'

    ! 解く行列をコピーする
    do i=1,nnode
        do j=1,nnode
            a(i,j)=syst(i,j)
            b(i,j)=sysu(i,j)
        end do
    end do

    ! eigrs(a,b,n,n1,ne,nv,eps,w,e,v,ier)
    !
    ! a(n1,n):  Ax=λBx の A                 (D)
    ! b(n1,n):  Ax=λBx の B                 (D)
    ! n:        行列の次数                   (I)
    ! n1:       配列の行数                   (I)
    ! ne:       求める固有値の数             (I)
    ! nv:       求める固有ベクトルの数       (I)
    ! eps:      収束判定基準                 (D)
    ! w(n1,7):  作業用領域                   (D)
    ! e(ne):    求められた固有値が入る       (D)
    ! v(n1,nv): 求められた固有ベクトルが入る (D)
    ! ier:      リターンコード               (I)
    !
    nne=nnode
    nnv=nne
    call eigab(a,b,nnode,nn1,-nne,nnv,1.0d-16,w,e,v,ier)

    ! 結果を変数に格納する
    do i=1,nne
        eigval(i)=e(i)
    end do

    do i=1,nnv
        do j=1,nnode
            eigvec(i,j)=v(j,i)
        end do
    end do

    return
end subroutine

!----------------------------------------------------------------------
! メッシュデータの確認
!----------------------------------------------------------------------
subroutine testmesh
   use mod_fem
   implicit none

   integer :: n(3),i,j
   real(8) :: x(3),y(3)

   write(*,*) '**** TESTMESH ****'

   open(10,file='mesh.dat')
   do i=1,nelem
      ! グローバル節点番号の取得
      do j=1,3
         n(j)=elem(i).node(j)
      end do

      ! x,y 座標の取得
      do j=1,3
         x(j)=gn2x(n(j))
         y(j)=gn2y(n(j))
      end do

      do j=1,3
         write(10,*) x(j),y(j)
      end do
      write(10,*) x(1),y(1) ! ポリゴンを閉じるため
      write(10,"()")        ! 改行
      write(10,"()")        ! 改行
   end do
   close(10)

   return
end subroutine

!----------------------------------------------------------------------
! 固有値をファイルに出力
!----------------------------------------------------------------------
subroutine outeig
    use mod_consts
    use mod_fem
    implicit none

    integer :: i,j,n

    write(*,*) '**** OUTEIG ****'

    !-------- 固有値の出力 --------
    ! 固有値はカットオフ波数の二乗(kc^2)となっている
    open(10,file='eigval.dat')
    do i=1,nnode
        write(10,*) i,eigval(i)
    end do
    close(10)

    !-------- カットオフ周波数 (GHz) の出力 --------
    open(10,file='fcu.dat')
    do i=1,nnode
        write(10,"(i3,f6.2)") &
           i,(velocity_of_light/(2.0d0*pi))*dsqrt(dabs(eigval(i)))/1.0d9
    end do
    close(10)

    return
end subroutine

!----------------------------------------------------------------------
! モード関数（ポテンシャルデータ）をファイルに出力
!----------------------------------------------------------------------
subroutine outpot(mode)
    use mod_fem
    implicit none

    integer mode

    integer :: n(3),i,j
    real(8) :: x(3),y(3)

    write(*,*) '**** OUTPOT ****'

    open(10,file='pot.dat')
    do i=1,nelem
        ! グローバル節点番号の取得
        do j=1,3
            n(j)=elem(i).node(j)
        end do

        ! x,y 座標の取得
        do j=1,3
            x(j)=gn2x(n(j))
            y(j)=gn2y(n(j))
        end do

        do j=1,3
            write(10,*) x(j),y(j),eigvec(mode,n(j))
        end do
        write(10,*) x(1),y(1),eigvec(mode,n(1)) ! ポリゴンを閉じるため
        write(10,"()")                          ! 改行
        write(10,"()")                          ! 改行
    end do
    close(10)

    return
end subroutine

!----------------------------------------------------------------------
! モード関数（フィールドデータ）をファイルに出力
!----------------------------------------------------------------------
subroutine outfield(mode)
    use mod_fem
    implicit none

    integer mode

    integer :: i,n1,n2,n3
    real(8) :: x1,y1,z1,x2,y2,z2,x3,y3,z3, &
               x(melem),y(melem),xc(melem),yc(melem), &
               nvecx,nvecy,nvecz,absmax

    write(*,*) '**** OUTFIELD ****'

    open(10,file='field.dat')
    absmax=0.0d0
    do i=1,nelem
        ! グローバル節点番号の取得
        n1=elem(i).node(1)
        n2=elem(i).node(2)
        n3=elem(i).node(3)

        ! x,y,z 座標の取得
        x1=gn2x(n1)
        y1=gn2y(n1)
        z1=eigvec(mode,n1)

        x2=gn2x(n2)
        y2=gn2y(n2)
        z2=eigvec(mode,n2)

        x3=gn2x(n3)
        y3=gn2y(n3)
        z3=eigvec(mode,n3)

        xc(i)=(x1+x2+x3)/3.0d0
        yc(i)=(y1+y2+y3)/3.0d0

        ! 3点を結ぶ法線ベクトルとその大きさ
        nvecx=-y2*z1+y3*z1+y1*z2-y3*z2-y1*z3+y2*z3
        nvecy= x2*z1-x3*z1-x1*z2+x3*z2+x1*z3-x2*z3
        nvecz=-x2*y1+x3*y1+x1*y2-x3*y2-x1*y3+x2*y3

        x(i)=nvecx
        y(i)=nvecy

        if(dsqrt(x(i)**2+y(i)**2).gt.absmax) then
            absmax=dsqrt(x(i)**2+y(i)**2)
        end if
    end do

    do i=1,nelem
        write(10,"(4e23.15)") &
                 xc(i),yc(i), &
                 xc(i)+y(i)/absmax*(sizey/ny),yc(i)+x(i)/absmax*(sizey/ny)
        write(10,"()")          ! 改行
        write(10,"()")          ! 改行
    end do
    close(10)

    return
end subroutine

!
! End of file
!