!**********************************************************************
! 有限要素法のテストルーチン
!**********************************************************************

!----------------------------------------------------------------------
! メッシュデータの確認 (2D 表示用)
!----------------------------------------------------------------------
subroutine testmesh
   use mod_fem
   implicit none

   integer :: nn(4),i,j,k
   real(8) :: x(4),y(4)

!   write(*,*) '**** TESTMESH-2D ****'
   write(1,*) '**** TESTMESH-2D ****'

   open(10,file='meshtest.dat')
   do i=1,nelem
      k=elem(i).npoly

      ! グローバル節点番号の取得
      do j=1,k
         nn(j)=elem(i).node(j)
      end do

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

      do j=1,k
         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

!----------------------------------------------------------------------
! 三角要素内の中心の電界を返す
!
! [Input]
! mode: モード番号
! ie: 要素番号（三角要素に限る）
!
! [Output]
! vx: 三角要素内の中心における電界の x 成分 Ex
! vy: 三角要素内の中心における電界の y 成分 Ey
!----------------------------------------------------------------------
subroutine field3c(mode,ie,vx,vy)
   use mod_fem
   implicit none

   integer :: mode,ie
   real(8) :: vx,vy

   integer :: edge(3),node1,node2,node3,i
   real(8) :: x1,x2,x3,y1,y2,y3,l(3),ll(3),a(3),b(3),c(3),area,&
              wx(3),wy(3),xc,yc
   real(8) :: dist

   ! 三角要素でないならエラー
   if(elem(ie).npoly.ne.3) then
      write(*,*) 'Sub(field3):'
      write(*,*) '   This element number do not represent triangle.'
      stop
   end if

   ! グローバルエッジへの変換
   edge(1)=elem(ie).edge(1)
   edge(2)=elem(ie).edge(2)
   edge(3)=elem(ie).edge(3)
   ! グローバルノードへの変換
   node1=elem(ie).node(1)
   node2=elem(ie).node(2)
   node3=elem(ie).node(3)
   ! 頂点座標の取得
   x1=gn2x(node1)
   x2=gn2x(node2)
   x3=gn2x(node3)

   y1=gn2y(node1)
   y2=gn2y(node2)
   y3=gn2y(node3)
   ! 中心座標の計算
   xc=(x1+x2+x3)/3.0d0
   yc=(y1+y2+y3)/3.0d0

   ! エッジの長さの計算
   l(1)=dist(x1,y1,x2,y2)
   l(2)=dist(x2,y2,x3,y3)
   l(3)=dist(x3,y3,x1,y1)
   ! パラメータ a_i,b_i,c_i
   a(1)=x2*y3-x3*y2
   a(2)=x3*y1-x1*y3
   a(3)=x1*y2-x2*y1
   b(1)=y2-y3
   b(2)=y3-y1
   b(3)=y1-y2
   c(1)=x3-x2
   c(2)=x1-x3
   c(3)=x2-x1
   ! 三角要素の面積
   area=(x1*b(1)+x2*b(2)+x3*b(3))/2.0d0
   ! 面積座標 L_i
   ll(1)=1.0d0/(2.0d0*area)*( a(1)+b(1)*xc+c(1)*yc )
   ll(2)=1.0d0/(2.0d0*area)*( a(2)+b(2)*xc+c(2)*yc )
   ll(3)=1.0d0/(2.0d0*area)*( a(3)+b(3)*xc+c(3)*yc )

   ! 基底ベクトル関数の計算
   wx(1)=l(1)/(2.0d0*area)*( ll(1)*b(2)-ll(2)*b(1) )
   wy(1)=l(1)/(2.0d0*area)*( ll(1)*c(2)-ll(2)*c(1) )

   wx(2)=l(2)/(2.0d0*area)*( ll(2)*b(3)-ll(3)*b(2) )
   wy(2)=l(2)/(2.0d0*area)*( ll(2)*c(3)-ll(3)*c(2) )

   wx(3)=l(3)/(2.0d0*area)*( ll(3)*b(1)-ll(1)*b(3) )
   wy(3)=l(3)/(2.0d0*area)*( ll(3)*c(1)-ll(1)*c(3) )

   ! 値を返す
   vx=0; vy=0
   do i=1,3
      vx=vx+eigvec(mode,edge(i))*( elem(ie).orient(i) )*wx(i)
      vy=vy+eigvec(mode,edge(i))*( elem(ie).orient(i) )*wy(i)
   end do

   return
end subroutine

!----------------------------------------------------------------------
! フィールドデータの確認
!
! [Input]
! mode: モード番号
!----------------------------------------------------------------------
subroutine testfield(mode)
   use mod_consts
   use mod_fem
   implicit none

   integer :: mode

   real(8) :: fieldxdat(melem),fieldydat(melem)
   integer :: n(4),i,j
   real(8) :: x(4),y(4),xc,yc,ampmax,amp,vabs,vx,vy,vvx,vvy,cnorm
   real(8) :: t=-90.0d0*(pi/180.0d0)
   vabs(vvx,vvy)=dsqrt(vvx**2+vvy**2)

!   write(*,*) '**** TESTFIELD ****'
   write(1,*) '**** TESTFIELD ****'

   ! フィールドデータを格納する
   do i=1,nelem
      if(elem(i).npoly.eq.3) then
         call field3c(mode,i,vx,vy)
      else
         write(*,*) "Sub(testfield): rectangular element is not supported!"
         stop
      end if
      fieldxdat(i)=vx
      fieldydat(i)=vy
      ! t [rad] だけ回転させる
      !fieldxdat(i)=dcos(t)*vx-dsin(t)*vy
      !fieldydat(i)=dsin(t)*vx+dcos(t)*vy
   end do

   ! フィールドの大きさの最大値を求める
   ampmax=0.0d0
   do i=1,nelem
      amp=vabs(fieldxdat(i),fieldydat(i))
      if(amp.gt.ampmax) then
         ampmax=amp
      end if
   end do
   ! フィールドデータを規格化する
   cnorm=min(sizex/nx,sizey/ny)/ampmax
   do i=1,nelem
      fieldxdat(i)=fieldxdat(i)*cnorm
      fieldydat(i)=fieldydat(i)*cnorm
   end do

   ! ファイルに書き込む
   open(10,file='fieldtest.dat')
   do i=1,nelem
      ! グローバル節点番号の取得
      do j=1,elem(i).npoly
         n(j)=elem(i).node(j)
      end do

      ! x,y 座標の取得
      do j=1,elem(i).npoly
         x(j)=gn2x(n(j))
         y(j)=gn2y(n(j))
      end do

      ! 要素の中心座標
      xc=0; yc=0;
      do j=1,elem(i).npoly
         xc=xc+x(j)
         yc=yc+y(j)
      end do
      xc=xc/elem(i).npoly
      yc=yc/elem(i).npoly

      write(10,"(4e23.15)") xc,yc,xc+fieldxdat(i),yc+fieldydat(i)
      write(10,"()")       ! 改行
      write(10,"()")       ! 改行
   end do

   close(10)

   return
end subroutine

!----------------------------------------------------------------------
! 要素番号の確認
!----------------------------------------------------------------------
subroutine testelemnum
   use mod_consts
   use mod_fem
   implicit none

   integer :: i,j,n(4)
   real(8) :: x(4),y(4),xc,yc

!   write(*,*) '**** TEST ELEM NUM ****'
   write(1,*) '**** TEST ELEM NUM ****'

   ! ファイルに書き込む
   open(10,file='elemnumtest.dat')
   do i=1,nelem
      ! グローバル節点番号の取得
      do j=1,elem(i).npoly
         n(j)=elem(i).node(j)
      end do

      ! x,y 座標の取得
      do j=1,elem(i).npoly
         x(j)=gn2x(n(j))
         y(j)=gn2y(n(j))
      end do

      ! 要素の中心座標
      xc=0; yc=0;
      do j=1,elem(i).npoly
         xc=xc+x(j)
         yc=yc+y(j)
      end do
      xc=xc/elem(i).npoly
      yc=yc/elem(i).npoly

      write(10,"(i,2f)") i,xc,yc
      write(10,"()")          ! 改行
      write(10,"()")          ! 改行
   end do

   close(10)

   return
end subroutine

!----------------------------------------------------------------------
! エッジ番号の確認
! mesh_orientation サブルーチン内でグローバルエッジ情報が作成
! されている必要がある
!----------------------------------------------------------------------
subroutine testedgenum
   use mod_consts
   use mod_fem
   implicit none

   integer :: i,j
   real(8) :: xc,yc

!   write(*,*) '**** TEST EDGE NUM ****'
   write(1,*) '**** TEST EDGE NUM ****'

   ! ファイルに書き込む
   open(10,file='edgenumtest.dat')
   do i=1,nedge
      ! エッジの中心座標
      xc=0.0d0; yc=0.0d0;
      do j=1,2
         xc=xc+gn2x(gedge(i).node(j))
         yc=yc+gn2y(gedge(i).node(j))
      end do
      xc=xc/2.0d0
      yc=yc/2.0d0

      write(10,"(i,2f)") i,xc,yc
      write(10,"()")          ! 改行
      write(10,"()")          ! 改行
   end do

   close(10)
   return
end subroutine

!----------------------------------------------------------------------
! ノード番号の確認
!----------------------------------------------------------------------
subroutine testnodenum
   use mod_consts
   use mod_fem
   implicit none

   integer :: i,j,n,npoly

!   write(*,*) '**** TEST NODE NUM ****'
   write(1,*) '**** TEST NODE NUM ****'

   ! ファイルに書き込む
   open(10,file='nodenumtest.dat')
   do i=1,nnode
      write(10,"(i,2f)") i,gn2x(i),gn2y(i)
      write(10,"()")          ! 改行
      write(10,"()")          ! 改行
   end do

   close(10)

   return
end subroutine

!----------------------------------------------------------------------
! 形状の確認
!
!   4                        3
!    ************************
!    ************************
!    ************************
!   1                        2
!
!----------------------------------------------------------------------
subroutine testshape
   use mod_fem
   implicit none

   integer,parameter :: ncorner=4
   integer :: corner(ncorner),i
   real(8) :: x(ncorner),y(ncorner)

!   write(*,*) '**** TEST SHAPE ****'
   write(1,*) '**** TEST SHAPE ****'

   corner(1)=elem(1).node(1)
   corner(2)=elem(nelem-2*ny+1).node(2)
   corner(3)=elem(nelem).node(3)
   corner(4)=elem(2*ny).node(1)

   open(10,file="testshape.dat")
   ! 導波管断面の描画
   ! x,y 座標の取得
   do i=1,ncorner
      x(i)=gn2x(corner(i))
      y(i)=gn2y(corner(i))
   end do

   do i=1,ncorner
      write(10,*) x(i),y(i)
   end do
   write(10,"(2(f))") x(1),y(1) ! ポリゴンを閉じるため
   write(10,"()")        ! 改行
   write(10,"()")        ! 改行

   return
end subroutine

!
! End of File
!