! gkaiten.f90やghanten.f90の要素分割を細かくしようとして
! 構造解析関係のプログラム置き場にあるccxs6s3.f90を使うと、
! なぜだか節点が重なってしまうところができるので、
! ccxs6s3.f90を真似て作り直してみる。
! つまり、S6(CCX)の要素数を4倍に増やすプログラム。

! 今は、螺旋折りプログラムに対応する用に、余計なものも書いておいてある。

!実行例：
! 元のinpファイルの
!   "*NODE, NSET=~" を消して節点数に、
!   "*ELEMENT, TYPE=S6, ELSET=~" を消して要素数に
! 書き換えて保存し、
! gfortran -o s6x4 s6x4.f90
! ./s6x4<old.inp>new.inp
! 材料定数・境界条件・荷重条件・出力の設定等もプログラムの中で
! 書き換えてからコンパイル・実行した方が楽。

implicit none
real :: x0(10000),y0(10000),z0(10000)
real :: x(30000),y(30000),z(30000)
integer :: mesh0(10000,6),mesh(30000,6),narabi(10000,10000)
integer i,j,k,l,nset0,nset,nyou0,nyou,kazu
integer nsaika

!#--螺旋折り---
   !#--入力
   real    :: r=32.85d-3      ! m   :半径
   real    :: h=120.d-3       ! m   :高さ
   real    :: t=2.d-4         ! m   :厚さ
   integer :: nc=5            !     :円周方向の分割数
   integer :: nh=5            !     :高さ方向の分割数
   real    :: p=1.d-3        ! MN  :荷重
   integer :: step=1          !     :計算/出力ステップ数
   !---
   real :: pi=2.d0*asin(1.d0) ! rad :円周率
   real dh,th,phi,hen1,hen2,hen3,alp,err
   dh=h/nh         ! m   :一段の高さ 
   phi=2.d0*pi/nc  ! rad :一辺の中心角
   !---一段ごとのずれる中心角(th)を計算
   do i=0,10000000!分割は適当
      th=pi*i/10000000!↑の分割に合わせる
      hen1=2.d0*r*sin(phi/2)
      hen2=sqrt((2.d0*r*sin(th/2))**2+dh**2)
      hen3=sqrt((2.d0*r*sin(th/2+phi/2))**2+dh**2)
      alp=acos((hen2**2+hen3**2-hen1**2)/(2.d0*hen2*hen3))
      err=alp-pi/nc
      if (err<1.d-12) exit!許容誤差は適当
   enddo
   print'(A,f12.8,A,i10)',"**th(deg)=",th*180.d0/pi,"  ，反復回数: ",i
   print'(A,f12.8,A,f12.8)',"**alpha(deg)=",alp*180.d0/pi," ≒  pi/m(deg)=",180.d0/nc
   !inpファイルに角度を出力
!#--螺旋折り---

narabi=-1!隣り合う節点の中間の節点番号の初期値をすべて-1にしておく
!---元inpファイルからデータ取得
read(*,*) nset0!元inpの節点数を読み取る
do i=1,nset0
   read(*,*) kazu,x0(i),y0(i),z0(i)!各節点の座標値の読み取り
   x(i)=x0(i)
   y(i)=y0(i)
   z(i)=z0(i)
enddo
read(*,*) nyou0!元inpの要素数を読み取る
do i=1,nyou0
   read(*,*) kazu,mesh0(i,1),mesh0(i,2),mesh0(i,3),&!各要素の節点番号の読み取り
             &mesh0(i,4),mesh0(i,5),mesh0(i,6)
enddo
!---新しく節点をふる
nset=nset0
do i=1,nyou0
   !---中間節点1
   if (narabi(mesh0(i,1),mesh0(i,4))==-1) then!既に節点をふってあるところは避ける
      nset=nset+1
      narabi(mesh0(i,1),mesh0(i,4))=nset!narabiに中間節点の番号を与える
      narabi(mesh0(i,4),mesh0(i,1))=nset!narabiに中間節点の番号を与える
      x(nset)=(x0(mesh0(i,1))+x0(mesh0(i,4)))/2
      y(nset)=(y0(mesh0(i,1))+y0(mesh0(i,4)))/2
      z(nset)=(z0(mesh0(i,1))+z0(mesh0(i,4)))/2
   endif
   !---中間節点2
   if (narabi(mesh0(i,4),mesh0(i,6))==-1) then
      nset=nset+1
      narabi(mesh0(i,4),mesh0(i,6))=nset
      narabi(mesh0(i,6),mesh0(i,4))=nset
      x(nset)=(x0(mesh0(i,4))+x0(mesh0(i,6)))/2
      y(nset)=(y0(mesh0(i,4))+y0(mesh0(i,6)))/2
      z(nset)=(z0(mesh0(i,4))+z0(mesh0(i,6)))/2
   endif
   !---中間節点3
   if (narabi(mesh0(i,6),mesh0(i,1))==-1) then
      nset=nset+1
      narabi(mesh0(i,6),mesh0(i,1))=nset
      narabi(mesh0(i,1),mesh0(i,6))=nset
      x(nset)=(x0(mesh0(i,6))+x0(mesh0(i,1)))/2
      y(nset)=(y0(mesh0(i,6))+y0(mesh0(i,1)))/2
      z(nset)=(z0(mesh0(i,6))+z0(mesh0(i,1)))/2
   endif
   !---中間節点4
   if (narabi(mesh0(i,4),mesh0(i,2))==-1) then
      nset=nset+1
      narabi(mesh0(i,4),mesh0(i,2))=nset
      narabi(mesh0(i,2),mesh0(i,4))=nset
      x(nset)=(x0(mesh0(i,4))+x0(mesh0(i,2)))/2
      y(nset)=(y0(mesh0(i,4))+y0(mesh0(i,2)))/2
      z(nset)=(z0(mesh0(i,4))+z0(mesh0(i,2)))/2
   endif
   !---中間節点5
   if (narabi(mesh0(i,2),mesh0(i,6))==-1) then
      nset=nset+1
      narabi(mesh0(i,2),mesh0(i,6))=nset
      narabi(mesh0(i,6),mesh0(i,2))=nset
      x(nset)=(x0(mesh0(i,2))+x0(mesh0(i,6)))/2
      y(nset)=(y0(mesh0(i,2))+y0(mesh0(i,6)))/2
      z(nset)=(z0(mesh0(i,2))+z0(mesh0(i,6)))/2
   endif
   !---中間節点6
   if (narabi(mesh0(i,2),mesh0(i,5))==-1) then
      nset=nset+1
      narabi(mesh0(i,2),mesh0(i,5))=nset
      narabi(mesh0(i,5),mesh0(i,2))=nset
      x(nset)=(x0(mesh0(i,2))+x0(mesh0(i,5)))/2
      y(nset)=(y0(mesh0(i,2))+y0(mesh0(i,5)))/2
      z(nset)=(z0(mesh0(i,2))+z0(mesh0(i,5)))/2
   endif
   !---中間節点7
   if (narabi(mesh0(i,5),mesh0(i,6))==-1) then
      nset=nset+1
      narabi(mesh0(i,5),mesh0(i,6))=nset
      narabi(mesh0(i,6),mesh0(i,5))=nset
      x(nset)=(x0(mesh0(i,5))+x0(mesh0(i,6)))/2
      y(nset)=(y0(mesh0(i,5))+y0(mesh0(i,6)))/2
      z(nset)=(z0(mesh0(i,5))+z0(mesh0(i,6)))/2
   endif
   !---中間節点8
   if (narabi(mesh0(i,5),mesh0(i,3))==-1) then
      nset=nset+1
      narabi(mesh0(i,5),mesh0(i,3))=nset
      narabi(mesh0(i,3),mesh0(i,5))=nset
      x(nset)=(x0(mesh0(i,5))+x0(mesh0(i,3)))/2
      y(nset)=(y0(mesh0(i,5))+y0(mesh0(i,3)))/2
      z(nset)=(z0(mesh0(i,5))+z0(mesh0(i,3)))/2
   endif
   !---中間節点9
   if (narabi(mesh0(i,3),mesh0(i,6))==-1) then
      nset=nset+1
      narabi(mesh0(i,3),mesh0(i,6))=nset
      narabi(mesh0(i,6),mesh0(i,3))=nset
      x(nset)=(x0(mesh0(i,3))+x0(mesh0(i,6)))/2
      y(nset)=(y0(mesh0(i,3))+y0(mesh0(i,6)))/2
      z(nset)=(z0(mesh0(i,3))+z0(mesh0(i,6)))/2
   endif
enddo
!---新しくふった節点座標を出力する
print'("*NODE, NSET=Nall")'
do i=1,nset
   print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
   &,i,',',x(i),',',y(i),',',z(i)
enddo
!---新しく要素を作り直す
nyou=0
do i=1,nyou0
   !---要素1
   nyou=nyou+1
   mesh(nyou,1)=mesh0(i,1)
   mesh(nyou,2)=mesh0(i,4)
   mesh(nyou,3)=mesh0(i,6)
   mesh(nyou,4)=narabi(mesh0(i,1),mesh0(i,4))
   mesh(nyou,5)=narabi(mesh0(i,4),mesh0(i,6))
   mesh(nyou,6)=narabi(mesh0(i,6),mesh0(i,1))
   !---要素2
   nyou=nyou+1
   mesh(nyou,1)=mesh0(i,6)
   mesh(nyou,2)=mesh0(i,4)
   mesh(nyou,3)=mesh0(i,2)
   mesh(nyou,4)=narabi(mesh0(i,6),mesh0(i,4))
   mesh(nyou,5)=narabi(mesh0(i,4),mesh0(i,2))
   mesh(nyou,6)=narabi(mesh0(i,2),mesh0(i,6))
   !---要素3
   nyou=nyou+1
   mesh(nyou,1)=mesh0(i,2)
   mesh(nyou,2)=mesh0(i,5)
   mesh(nyou,3)=mesh0(i,6)
   mesh(nyou,4)=narabi(mesh0(i,2),mesh0(i,5))
   mesh(nyou,5)=narabi(mesh0(i,5),mesh0(i,6))
   mesh(nyou,6)=narabi(mesh0(i,6),mesh0(i,2))
   !---要素4
   nyou=nyou+1
   mesh(nyou,1)=mesh0(i,6)
   mesh(nyou,2)=mesh0(i,5)
   mesh(nyou,3)=mesh0(i,3)
   mesh(nyou,4)=narabi(mesh0(i,6),mesh0(i,5))
   mesh(nyou,5)=narabi(mesh0(i,5),mesh0(i,3))
   mesh(nyou,6)=narabi(mesh0(i,3),mesh0(i,6))
enddo
!---新しく作った要素の出力
print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
do i=1,nyou
   print'(i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8)',&
       &i,',',&
       &mesh(i,1),',',&
       &mesh(i,2),',',&
       &mesh(i,3),',',&
       &mesh(i,4),',',&
       &mesh(i,5),',',&
       &mesh(i,6)
end do

!#--螺旋折り---
   !---固定端の節点群を出力
   print'("*NSET,NSET=kotei")'
   do i=1,nset
      if (z(i)>-1.d-6.and.z(i)<1.d-6) then
         print'(i8,a)',i
      endif
   enddo
   !---載荷点の節点群を出力
   nsaika=0
   print'("*NSET,NSET=saika")'
   do i=1,nset
      if (z(i)>h-1.d-6.and.z(i)<h+1.d-6) then
         print'(i8,a)',i
         nsaika=nsaika+1
      endif
   enddo
   !print'("*RIGID BODY,NSET=saika")'!載荷部分は平らに変形させる(強制変位では使えない)
   !---段ごとのz変位を均等にする
   !print'("*EQUATION")'
   !do j=1,nh-1
   !   do i=1,nset
   !      if (z(i)>j*dh-1.d-6.and.z(i)<j*dh+1.d-6) then
   !         print'("2")'!項が２つという意味
   !         print'(i7,A,f10.7,a,i7,A,f10.7)'&
   !               &,i,',3,',real(nh)/real(j),',',上の節点のどこか,',3,',-1.d0!ちょっと面倒臭い
   !      endif
   !   enddo
   !enddo
   !#--境界条件
   print'("*BOUNDARY")'
   print'("kotei,1,3")'
   !#---材料定数と厚さ
   print'("*MATERIAL,NAME=EL1")'
   print'("*ELASTIC")'
   print'("69d3, 0.3d0")'!アルミ缶69GPa
   print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL1")'
   print'(1pe13.5,A,A)',t,',,','EL1'
   !---1ステップごとの計算/出力
   do i=1,step
      !#--載荷条件
      print'("*STEP,NLGEOM,INC=50000")'!幾何学非線形の場合
      !print'("*STEP")'!線形の場合
      print'("*STATIC,SOLVER=SPOOLES")'
      !print'("0.001,1")'!幾何学非線形のなんとか増分
      !print'("*BUCKLE")'!座屈計算の場合
      !print'("10,0.01")'
      print'("*CLOAD")'!載荷
      print'(A,a,i6,a,1pe13.5)','saika',',',3,',',-p/nsaika
      !print'("*BOUNDARY")'!強制変位
      !print'(A,a,i6,A,1pe13.5)','saika',',',3,',,',(-h+t*nh)*i/step!厚みを考慮して高さ0まで畳む
      !---出力の設定
      print'("*NODE PRINT,NSET=saika")'
      print'("U,RF")'
      !print'("*EL PRINT,ELSET=Eall")'
      !print'("S")'
      print'("*NODE FILE")'
      print'("U")'
      print'("*EL FILE")'
      print'("S")'
      print'("*END STEP")'
   enddo
!#--螺旋折り---

end
