!frdファイル(CCX:C3D10)をvtkファイル(paraview)に変換する。

!使い方:frdファイルを開く。frdファイルは'-'の前に空欄がないので、
!       :%s/-/\ -/g でEnter、:%s/E\ -/E-/g でEnterで、保存。
!       計算ステップ数と材料数と入出力ファイル名を書き換えて実行。
!       ただし，frdファイルにひずみ(*EL PRINT の E)を出力していない場合は，
!       ひずみに関する行をコメントアウトしておく．
!       $ paraview hoge..vtk& で全ステップのファイルを開ける。

implicit none
!#--手入力
integer :: nsteps=1               !計算ステップ数
integer :: nmaterials=2           !材料数
!---
integer                nnodes     !節点数
integer                nvolumes   !要素数
real(8),allocatable :: x(:)       !元座標
real(8),allocatable :: y(:)
real(8),allocatable :: z(:)
integer,allocatable :: youso(:,:) !要素番号とその節点番号
real(8) u,v,w !節点変位
real(8),allocatable :: sxx(:)     !応力(true Cauchy stress)
real(8),allocatable :: syy(:)
real(8),allocatable :: szz(:)
real(8),allocatable :: sxy(:)
real(8),allocatable :: syz(:)
real(8),allocatable :: szx(:)
real(8),allocatable :: exx(:)     !ひずみ(total Lagrangian strain)
real(8),allocatable :: eyy(:)
real(8),allocatable :: ezz(:)
real(8),allocatable :: exy(:)
real(8),allocatable :: eyz(:)
real(8),allocatable :: ezx(:)
real(8) :: vmis    !von Mises応力
real(8) :: a(3,3)  !応力の行列化(主応力の計算に使用)
real(8) :: vc(3,3)  !固有ベクトル
real(8) time  !時間
integer itime !整数(時間*10^3)
integer t,i,j,k,icut,n
character :: ccut*20,filename*30,m*10

!#--入力ファイル名
open(11,file='c3d10.frd',status='old')!読み取り専用

!---節点変位の読み込み
do i=1,8+nmaterials
   read(11,'()') 
enddo
read(11,*) ccut,nnodes
allocate(x(nnodes),y(nnodes),z(nnodes))!配列の確保
do i=1,nnodes
   read(11,*) icut,icut,x(i),y(i),z(i)
enddo

!---要素情報の読み込み
read(11,'()')
read(11,*) ccut,nvolumes
allocate(youso(nvolumes,10))!配列の確保
do i=1,nvolumes
   read(11,'()')
   read(11,*) icut,youso(i,1),youso(i,2),youso(i,3),youso(i,4),youso(i,5),&
       &youso(i,6),youso(i,7),youso(i,8),youso(i,9),youso(i,10)
enddo

!#--出力ファイル名
do t=1,nsteps
   time=real(t)/nsteps
   itime=int(time*1.d3)
   write(m,'(i10)') itime
   filename='post/c3d10_'//trim(adjustl(m))//'.vtk'
   open(12,file=(filename),status='replace')!新規置き換え
   write(12,'(A)') '# vtk DataFile Version 2.0'
   write(12,'(A)') 'frd2vtk'!ヘッダー:何でもいいから書かないといけない
   write(12,'(A)') 'ASCII'
   write(12,'(A)') 'DATASET UNSTRUCTURED_GRID'
!---節点座標の書き込み
   write(12,'(A,i10,A)') 'POINTS ',nnodes,' float'
   do i=1,nnodes
      write(12,'(1pe17.9,1pe17.9,1pe17.9)') x(i),y(i),z(i)
   enddo
!---要素の並び順の書き込み
   write(12,'(A,i10,i15)') 'CELLS ',nvolumes,nvolumes*(10+1)
   do i=1,nvolumes
      write(12,'(i2,i8,i8,i8,i8,i8,i8,i8,i8,i8,i8)') 10,&
          &youso(i,1)-1,youso(i,2)-1,youso(i,3)-1,youso(i,4)-1,youso(i,5)-1,&
          &youso(i,6)-1,youso(i,7)-1,youso(i,8)-1,youso(i,9)-1,youso(i,10)-1
   enddo
!---要素の形状の書き込み
   write(12,'(A,i8)') 'CELL_TYPES ',nvolumes
   do i=1,nvolumes
      write(12,*) '24' !24はC3D10のCELL_TYPE
   enddo
enddo
deallocate(x,y,z,youso)!動的メモリの解放

do t=1,nsteps
!#--vtkファイルを再び開く
   time=real(t)/nsteps
   itime=int(time*1.d3)
   write(m,'(i10)') itime
   filename='post/c3d10_'//trim(adjustl(m))//'.vtk'!名前を上のvtkファイルに合わせる
   open(12,file=(filename),position='append')!途中から書き足す

!---節点変位の読み書き
   do i=1,8
      read(11,'()')
   enddo
   write(12,'(A,i10)') 'POINT_DATA ',nnodes
   write(12,'(A)') 'VECTORS DXDYDZ float'!DXDYDZは何でもいい
   do i=1,nnodes
      read(11,*) icut,icut,u,v,w
      write(12,'(1pe17.9,1pe17.9,1pe17.9)') u,v,w
   enddo
!---応力テンソルの読み書き
   allocate(sxx(nnodes),syy(nnodes),szz(nnodes))!配列の確保
   allocate(sxy(nnodes),syz(nnodes),szx(nnodes))
   do i=1,10
      read(11,'()')
   enddo
   do i=1,nnodes
      read(11,*) icut,icut,sxx(i),syy(i),szz(i),sxy(i),syz(i),szx(i)
   enddo
   write(12,'(A)') 'SCALARS SXX float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') sxx(i)
   enddo
   write(12,'(A)') 'SCALARS SYY float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') syy(i)
   enddo
   write(12,'(A)') 'SCALARS SZZ float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') szz(i)
   enddo
   write(12,'(A)') 'SCALARS SXY float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') sxy(i)
   enddo
   write(12,'(A)') 'SCALARS SYZ float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') syz(i)
   enddo
   write(12,'(A)') 'SCALARS SZX float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') szx(i)
   enddo
!---主応力の書き出し
   n=3
   write(12,'(A)') 'VECTORS PRIN1 float'
   do i=1,nnodes   
      a(1,:)=(/ sxx(i),sxy(i),szx(i) /)
      a(2,:)=(/ sxy(i),syy(i),syz(i) /)
      a(3,:)=(/ szx(i),syz(i),szz(i) /)
      call sort_eigen_jacobi(n,a,vc)
      write(12,'(1pe17.9,1pe17.9,1pe17.9)') vc(1,3)*a(3,3),vc(2,3)*a(3,3),vc(3,3)*a(3,3)
   enddo
   write(12,'(A)') 'VECTORS PRIN2 float'
   do i=1,nnodes   
      a(1,:)=(/ sxx(i),sxy(i),szx(i) /)
      a(2,:)=(/ sxy(i),syy(i),syz(i) /)
      a(3,:)=(/ szx(i),syz(i),szz(i) /)
      call sort_eigen_jacobi(n,a,vc)
      write(12,'(1pe17.9,1pe17.9,1pe17.9)') vc(1,2)*a(2,2),vc(2,2)*a(2,2),vc(3,2)*a(2,2)
   enddo
   write(12,'(A)') 'VECTORS PRIN3 float'
   do i=1,nnodes   
      a(1,:)=(/ sxx(i),sxy(i),szx(i) /)
      a(2,:)=(/ sxy(i),syy(i),syz(i) /)
      a(3,:)=(/ szx(i),syz(i),szz(i) /)
      call sort_eigen_jacobi(n,a,vc)
      write(12,'(1pe17.9,1pe17.9,1pe17.9)') vc(1,1)*a(1,1),vc(2,1)*a(1,1),vc(3,1)*a(1,1)
   enddo
!---von Mises応力の書き出し
   write(12,'(A)') 'SCALARS VMIS float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      vmis=sqrt(((sxx(i)-syy(i))**2.d0+(syy(i)-szz(i))**2.d0+(szz(i)-sxx(i))**2.d0&
          &+6.0d0*(sxy(i)**2.d0+syz(i)**2.d0+szx(i)**2.d0))/2.d0)
      write(12,'(1pe17.9)') vmis
   enddo
   deallocate(sxx,syy,szz,sxy,syz,szx)!動的メモリの解放
!#--ひずみの読み書き
   allocate(exx(nnodes),eyy(nnodes),ezz(nnodes))!配列の確保
   allocate(exy(nnodes),eyz(nnodes),ezx(nnodes))
   do i=1,10
      read(11,'()')
   enddo
   do i=1,nnodes
      read(11,*) icut,icut,exx(i),eyy(i),ezz(i),exy(i),eyz(i),ezx(i)
   enddo
   write(12,'(A)') 'SCALARS EXX float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') exx(i)
   enddo
   write(12,'(A)') 'SCALARS EYY float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') eyy(i)
   enddo
   write(12,'(A)') 'SCALARS EZZ float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') ezz(i)
   enddo
   write(12,'(A)') 'SCALARS EXY float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') exy(i)
   enddo
   write(12,'(A)') 'SCALARS EYZ float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') eyz(i)
   enddo
   write(12,'(A)') 'SCALARS EZX float'
   write(12,'(A)') 'LOOKUP_TABLE default'
   do i=1,nnodes
      write(12,'(1pe17.9)') ezx(i)
   enddo
   close(12)
   deallocate(exx,eyy,ezz,exy,eyz,ezx)!動的メモリの解放
enddo
close(11)
end

! 実対称行列の固有値をヤコビ法で求めた後，昇順に並び替えるサブルーチン
! ヤコビ法については，
!   戸川隼人『詳解 数値計算演習』(共立出版)p.163-169
! を参考にした．
subroutine sort_eigen_jacobi(n,a,v)
  implicit none
  integer,intent(in)    :: n      ! 行列の大きさ
  real(8),intent(inout) :: a(n,n) ! 対象とする実対称行列
  real(8),intent(out)   :: v(n,n) ! 固有ベクトル
  real(8) :: a1(n,n)              ! a_k+1
  real(8) :: v1(n,n)              ! v_k+1
  real(8) :: err=1.d-8            ! 許容誤差
  real(8) :: amax,th,c,s
  real(8) :: akari                ! 固有値を並び替えるときに使用
  real(8) :: vkari(n)             ! 固有ベクトルを並び替えるときに使用
  integer :: i,j,k,l

  a1=a
  ! 初期の固有ベクトルの単位行列化
  v=0.d0
  do i=1,n
     v(i,i)=1.d0
  enddo
  v1=v

  ! 反復計算のループ
  do
    ! 非対角行列の中で絶対値が最大なものを探す
    amax=0.d0
    do j=2,n
       do i=1,j-1
          if (abs(a(i,j))>amax) then
             amax=abs(a(i,j))
             k=i
             l=j
          endif
       enddo
    enddo
    if (amax**2<err) exit ! 収束判定(誤差はamax^2のオーダー)

    ! 行列P中のsin(th),cos(th)
    th=atan(2.d0*a(k,l)/(a(k,k)-a(l,l)))/2
    c=cos(th)
    s=sin(th)

    ! PAP^tの計算(効率的にk,l行とk,l列だけを計算する)
    a1(k,:)=a(k,:)*c+a(l,:)*s
    a1(l,:)=a(l,:)*c-a(k,:)*s
    a1(:,k)=a1(k,:)
    a1(:,l)=a1(l,:)
    a1(k,k)=a(k,k)*c**2+2.d0*a(k,l)*c*s+a(l,l)*s**2
    a1(l,l)=a(l,l)*c**2-2.d0*a(k,l)*c*s+a(k,k)*s**2
    a1(k,l)=a(k,l)*c**2+(a(l,l)-a(k,k))*c*s-a(k,l)*s**2
    a1(l,k)=a(k,l)*c**2+(a(l,l)-a(k,k))*c*s-a(k,l)*s**2
    a=a1 ! Aの更新

    ! VP^tの計算(効率的にi,j列だけを計算する)
    v1(:,k)=v(:,k)*c+v(:,l)*s
    v1(:,l)=v(:,l)*c-v(:,k)*s
    v=v1 ! Vの更新 

  enddo

  ! 昇順に並び替える
  do i=1,n-1
     do j=1,n-i
        if (a(j,j)>a(j+1,j+1)) then
           akari=a(j+1,j+1)
           a(j+1,j+1)=a(j,j)
           a(j,j)=akari
           vkari=v(:,j+1)
           v(:,j+1)=v(:,j)
           v(:,j)=vkari
        endif
     enddo
  enddo
end
