!PC鋼棒のねじれと摩擦(CCX)
! 円柱(B32)+円盤(C3D6+C3D8)*円盤(C3D6+C3D8) の接触
! PC鋼棒に引張力(p)と導入トルク(t)を与える。
! 円盤の大きさや摩擦係数は、何か根拠を持たせたい。

! datファイルの見方や出力の設定は、わからない点も多いが、
! KNOT
! tra  節点番号   u     v     w
! rot      ?    th_x  th_y  th_z
! edp      ?      ?     ?     ?

implicit none
!#--入力
real    :: d=13.d-3    ! m    : 直径
real    :: ell=1.825d0 ! m    : 長さ
real    :: t=230.d-6   ! MNm  : 導入トルク
real    :: p=70.d-6    ! MN   : 引張力
integer :: nstep=100   !      : 出力ステップ
integer :: nz=1000     !      : 軸方向分割数
real    :: r=65.d-3    ! m    : 摩擦用円盤の半径
integer :: nenban=4    ! 偶数 : 摩擦用円盤の高さ方向要素数
integer :: nr=8        !      : 中心方向分割数
integer :: nc=32       !      : 接線方向分割数
real    :: e1=206.d3   ! MPa  : PC鋼棒のヤング率
real    :: nu1=0.3d0   !      : PC鋼棒のポアソン比
real    :: e2=206.d3   ! MPa  : 円盤のヤング率(大きめ)
real    :: nu2=0.499d0 !      : 円盤のポアソン比
real    :: mu=0.15d4   !      : 円盤の摩擦係数
!---
real    :: dz,pi,dth,dr,ts
real    :: x(10000),y(10000),z(10000)
integer :: mesh0(100,3),mesh1(10000,3),mesh2(10000,6),mesh3(10000,8)
integer :: nset,nset1,nset2,nset3,nyou,nyou0,nyou1,nyou2,nyou3,nmen
integer :: i,j,k,l
!---
pi=2.d0*asin(1.d0) ! rad  : 円周率
dth=pi*2.d0/nc     ! rad  : 円周方向の要素角度
dr=r/nr        ! m    : 半径方向の要素長さ
dz=ell/nz          ! m    : 軸方向の要素長さ
nmen=nr*nc         !      : 円盤断面の節点数-1
!---節点座標
nset=0
do k=0,nz*2!B32
   nset=nset+1
   x(nset)=0.d0
   y(nset)=0.d0
   z(nset)=k*dz/2
end do
nset1=nset
do k=0,nenban!C3D6 & C3D8 1
   do j=1,nr
      do i=0,nc-1
         nset=nset+1
         x(nset)=dr*j*cos(dth*i)
         y(nset)=dr*j*sin(dth*i)
         z(nset)=dz/2*k
      end do
   end do
enddo
nset2=nset-nset1
do k=0,nenban!C3D6 & C3D8 1
   nset=nset+1
   x(nset)=0.d0
   y(nset)=0.d0
   z(nset)=dz/2*(nenban+k)
   do j=1,nr
      do i=0,nc-1
         nset=nset+1
         x(nset)=dr*j*cos(dth*i)
         y(nset)=dr*j*sin(dth*i)
         z(nset)=dz/2*(nenban+k)
      end do
   end do
enddo
nset3=nset-nset2-nset1
!---節点座標の出力
print'("*NODE, NSET=Nall")'
do i=1,nset
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
    &i,',',x(i),',',y(i),',',z(i)
end do
nyou1=0
do k=0,nz-1!B32
   nyou1=nyou1+1
   mesh1(nyou1,1)=k*2+1
   mesh1(nyou1,2)=mesh1(nyou1,1)+1
   mesh1(nyou1,3)=mesh1(nyou1,2)+1
end do
nyou0=0
do i=1,nc!B32 と C3D6 1 を一体化
   nyou0=nyou0+1
   mesh0(nyou0,1)=nenban/2+1
   mesh0(nyou0,2)=nset1+nmen*nenban/2+i
   mesh0(nyou0,3)=mesh0(nyou0,2)+nc
enddo
nyou2=0
nyou3=0
do k=0,nenban-1!C3D6 1
   do i=0,nc-2
   nyou2=nyou2+1
   mesh2(nyou2,1)=1+k
   mesh2(nyou2,2)=nset1+nmen*k+i+1
   mesh2(nyou2,3)=mesh2(nyou2,2)+1
   mesh2(nyou2,4)=mesh2(nyou2,1)+1
   mesh2(nyou2,5)=mesh2(nyou2,2)+nmen
   mesh2(nyou2,6)=mesh2(nyou2,3)+nmen
   end do
   nyou2=nyou2+1
   mesh2(nyou2,1)=1+k
   mesh2(nyou2,2)=nset1+nmen*k+nc
   mesh2(nyou2,3)=mesh2(nyou2,2)-nc+1
   mesh2(nyou2,4)=mesh2(nyou2,1)+1
   mesh2(nyou2,5)=mesh2(nyou2,2)+nmen
   mesh2(nyou2,6)=mesh2(nyou2,3)+nmen
   do j=0,nr-2!C3D8 1
      do i=1,nc-1
         nyou3=nyou3+1
         mesh3(nyou3,1)=nset1+nmen*k+nc*j+i
         mesh3(nyou3,2)=mesh3(nyou3,1)+nc
         mesh3(nyou3,3)=mesh3(nyou3,2)+1
         mesh3(nyou3,4)=mesh3(nyou3,3)-nc
         mesh3(nyou3,5)=mesh3(nyou3,1)+nmen
         mesh3(nyou3,6)=mesh3(nyou3,2)+nmen
         mesh3(nyou3,7)=mesh3(nyou3,3)+nmen
         mesh3(nyou3,8)=mesh3(nyou3,4)+nmen
      end do
   nyou3=nyou3+1
   mesh3(nyou3,1)=nset1+nmen*k+nc*j+i
   mesh3(nyou3,2)=mesh3(nyou3,1)+nc
   mesh3(nyou3,3)=mesh3(nyou3,2)+1-nc
   mesh3(nyou3,4)=mesh3(nyou3,3)-nc
   mesh3(nyou3,5)=mesh3(nyou3,1)+nmen
   mesh3(nyou3,6)=mesh3(nyou3,2)+nmen
   mesh3(nyou3,7)=mesh3(nyou3,3)+nmen
   mesh3(nyou3,8)=mesh3(nyou3,4)+nmen
   end do
end do
do k=0,nenban-1!C3D6 2
   do i=1,nc-1
      nyou2=nyou2+1
      mesh2(nyou2,1)=nset1+nset2+(nmen+1)*k+1
      mesh2(nyou2,2)=mesh2(nyou2,1)+i
      mesh2(nyou2,3)=mesh2(nyou2,2)+1
      mesh2(nyou2,4)=mesh2(nyou2,1)+(nmen+1)
      mesh2(nyou2,5)=mesh2(nyou2,2)+(nmen+1)
      mesh2(nyou2,6)=mesh2(nyou2,3)+(nmen+1)
   end do
   nyou2=nyou2+1
   mesh2(nyou2,1)=nset1+nset2+(nmen+1)*k+1
   mesh2(nyou2,2)=mesh2(nyou2,1)+i
   mesh2(nyou2,3)=mesh2(nyou2,2)-nc+1
   mesh2(nyou2,4)=mesh2(nyou2,1)+(nmen+1)
   mesh2(nyou2,5)=mesh2(nyou2,2)+(nmen+1)
   mesh2(nyou2,6)=mesh2(nyou2,3)+(nmen+1)
   do j=0,nr-2!C3D8 2
      do i=1,nc-1
         nyou3=nyou3+1
         mesh3(nyou3,1)=nset1+nset2+1+(nmen+1)*k+nc*j+i
         mesh3(nyou3,2)=mesh3(nyou3,1)+nc
         mesh3(nyou3,3)=mesh3(nyou3,2)+1
         mesh3(nyou3,4)=mesh3(nyou3,3)-nc
         mesh3(nyou3,5)=mesh3(nyou3,1)+(nmen+1)
         mesh3(nyou3,6)=mesh3(nyou3,2)+(nmen+1)
         mesh3(nyou3,7)=mesh3(nyou3,3)+(nmen+1)
         mesh3(nyou3,8)=mesh3(nyou3,4)+(nmen+1)
      end do
   nyou3=nyou3+1
   mesh3(nyou3,1)=nset1+nset2+1+(nmen+1)*k+nc*j+nc
   mesh3(nyou3,2)=mesh3(nyou3,1)+nc
   mesh3(nyou3,3)=mesh3(nyou3,2)+1-nc
   mesh3(nyou3,4)=mesh3(nyou3,3)-nc
   mesh3(nyou3,5)=mesh3(nyou3,1)+(nmen+1)
   mesh3(nyou3,6)=mesh3(nyou3,2)+(nmen+1)
   mesh3(nyou3,7)=mesh3(nyou3,3)+(nmen+1)
   mesh3(nyou3,8)=mesh3(nyou3,4)+(nmen+1)
   end do
end do
!---要素の出力
print'("*ELEMENT, TYPE=B32, ELSET=E1")'
do i=1,nyou1
print'(i6,a,i6,a,i6,a,i6)',&
    &i,',',&
    &mesh1(i,1),',',&
    &mesh1(i,2),',',&
    &mesh1(i,3)
end do
print'("*ELEMENT, TYPE=B32, ELSET=E0")'
do i=1,nyou0
print'(i6,a,i6,a,i6,a,i6)',&
    &nyou1+i,',',&
    &mesh0(i,1),',',&
    &mesh0(i,2),',',&
    &mesh0(i,3)
end do
print'("*ELEMENT, TYPE=C3D6, ELSET=E2")'
do i=1,nyou2
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &nyou0+nyou1+i,',',&
    &mesh2(i,1),',',&
    &mesh2(i,2),',',&
    &mesh2(i,3),',',&
    &mesh2(i,4),',',&
    &mesh2(i,5),',',&
    &mesh2(i,6)
end do
print'("*ELEMENT, TYPE=C3D8, ELSET=E3")'
do i=1,nyou3
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &nyou0+nyou1+nyou2+i,',',&
    &mesh3(i,1),',',&
    &mesh3(i,2),',',&
    &mesh3(i,3),',',&
    &mesh3(i,4),',',&
    &mesh3(i,5),',',&
    &mesh3(i,6),',',&
    &mesh3(i,7),',',&
    &mesh3(i,8)
end do
!---接触面の定義
print'(A)',"*SURFACE,NAME=SM"!マスター面群(ネジと一体化している円盤)
do i=1,nc!C3D6
   l=nyou0+nyou1+nyou2/2-nc+i
   print'(i6,A)',l,',S2'
enddo
l=nyou0+nyou1+nyou2+nyou3/2-((nr-1)*nc)
do j=0,nr-2!C3D8
   do i=1,nc
      l=l+1
      print'(i6,A)',l,',S2'
   enddo
enddo
print'(A)',"*SURFACE,NAME=SS"!スレーブ面群(宙に浮いている方の円盤)
do i=1,nc!C3D6
   l=nyou0+nyou1+nyou2/2+i
   print'(i6,A)',l,',S1'
enddo
l=nyou0+nyou1+nyou2+nyou3/2
do j=0,nr-2!C3D8
   do i=1,nc
      l=l+1
      print'(i6,A)',l,',S1'
   enddo
enddo
!---境界条件と荷重条件の節点群
print'("*NSET,NSET=saika")'
print'(i6)',nset1
print'("*NSET,NSET=ziku")'
print'(i6)',1
print'("*NSET,NSET=kotei")'
l=nset1+nset2+nset3-(nc*nr)
print'(i6)',l
do j=1,nr
   do i=1,nc
      l=l+1
      print'(i6)',l
   enddo
enddo
!#--接触と摩擦の設定
!print'(A)',"*CONTACT PAIR,INTERACTION=SI,TYPE=NODE TO SURFACE"!ccx_2.9
print'(A)',"*CONTACT PAIR,INTERACTION=SI"!ccx_2.6
print'(A)',"SS,SM"!スレーブ面,マスター面
print'(A)',"*SURFACE INTERACTION,NAME=SI"
print'("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=EXPONENTIAL")'
print'("1.e-4,1.e-3")'!接触と判定する距離(m),距離0の時の応力?
print'("*FRICTION")'
print'(1pe13.5,a,i6)',mu,",",1000!摩擦係数,謎パラメータ
!#--材料設定
print'("*MATERIAL,NAME=EL1")'!PC鋼棒
print'("*ELASTIC")'
print'(1pe13.5,a,f12.4)',e1,',',nu1
print'("*MATERIAL,NAME=EL2")'!円盤
print'("*ELASTIC")'
print'(1pe13.5,a,f12.4)',e2,',',nu2
print'("*BEAM SECTION,ELSET=E1,MATERIAL=EL1,SECTION=CIRC")'
print'(f10.6,a,f10.6)',d,',',d
print'(f10.6,a,f10.6,a,f10.6)',0.d0,',',1.d0,',',0.d0
print'("*BEAM SECTION,ELSET=E0,MATERIAL=EL2,SECTION=RECT")'
print'(f10.6,a,f10.6)',dz/2*nenban,',',dz/2*nenban
print'(f10.6,a,f10.6,a,f10.6)',0.d0,',',1.d0,',',0.d0
print'("*SOLID SECTION,ELSET=E2,MATERIAL=EL2")'
print'("*SOLID SECTION,ELSET=E3,MATERIAL=EL2")'
!#--境界条件
print'("*BOUNDARY")'
print'("kotei,1,3")'
print'("ziku,1,2")'
print'("ziku,4,5")'
print'("saika,1,2")'
print'("saika,4,5")'
!---出力ステップ
print'("*TIME POINTS,NAME=T1")'
do i=1,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo
!---載荷条件(引張)
print'("*STEP")'
print'("*STATIC")'
!print'("0.1,1.")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',p
!---datのオプション
print'("*NODE PRINT,NSET=saika")'
print'("U")'
print'("*NODE PRINT,NSET=ziku")'
print'("U")'
!print'("*EL PRINT,ELSET=E1")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
print'("S")'
print'("*END STEP")'
!---載荷条件(ねじり)
print'("*STEP,NLGEOM,INC=1000")'
print'("*STATIC")'
!print'("0.1,1.")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',6,',',t
!---datのオプション
print'("*NODE PRINT,NSET=saika,TIMEPOINTS=T1")'
print'("U")'
print'("*NODE PRINT,NSET=ziku,TIMEPOINTS=T1")'
print'("U")'
!print'("*EL PRINT,ELSET=E1,TIMEPOINTS=T1")'
!print'("S")'
print'("*NODE FILE,TIMEPOINTS=T1")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T1")'
print'("S")'
print'("*END STEP")'

end
