! gkaiten.f90繧�ghanten.f90縺ｮ隕∫ｴ�蛻�牡繧堤ｴｰ縺九￥縺励ｈ縺�→縺励※
! 讒矩�隗｣譫宣未菫ゅ�繝励Ο繧ｰ繝ｩ繝�鄂ｮ縺榊�ｴ縺ｫ縺ゅｋccxs6s3.f90繧剃ｽｿ縺�→縲�
! 縺ｪ縺懊□縺狗ｯ轤ｹ縺碁㍾縺ｪ縺｣縺ｦ縺励∪縺�→縺薙ｍ縺後〒縺阪ｋ縺ｮ縺ｧ縲�
! ccxs6s3.f90繧堤悄莨ｼ縺ｦ菴懊ｊ逶ｴ縺励※縺ｿ繧九�
! 縺､縺ｾ繧翫ヾ6(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縺ｫ縺励※縺翫￥
!---蜈ナnp繝輔ぃ繧､繝ｫ縺九ｉ繝��繧ｿ蜿門ｾ�
read(*,*) nset0!蜈ナnp縺ｮ遽轤ｹ謨ｰ繧定ｪｭ縺ｿ蜿悶ｋ
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!蜈ナnp縺ｮ隕∫ｴ�謨ｰ繧定ｪｭ縺ｿ蜿悶ｋ
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