!折り目を剛結にした反転螺旋折り円筒をCCX用に作ってみる
! 板要素(S6)で、随分と粗い
implicit none
!#--入力
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=1000.d-6                    ! MN  :荷重
integer :: step=10                       !     :計算/出力ステップ数
!---
real :: x(1000000),y(1000000),z(1000000) ! m   :節点の座標
integer :: mesh(1000000,6)               !     :板要素と立体要素の節点群
real :: pi=2.d0*asin(1.d0)               ! rad :円周率
integer i,j,k,nset,nset1,ntyo,ntei,nsya,nyou
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ファイルに角度を出力

!---頂点の節点座標
nset=0
do j=0,nh
   do i=0,nc-1
      nset=nset+1
      if (mod(j,2)==1) then
         x(nset)=r*cos(phi*i+th)
         y(nset)=r*sin(phi*i+th)
         z(nset)=j*dh
      else
         x(nset)=r*cos(phi*i)
         y(nset)=r*sin(phi*i)
         z(nset)=j*dh
      endif
   enddo
enddo
ntyo=nset !頂点の節点数
!---底辺の中間節点の座標
nset1=0
do j=0,nh
   do i=0,nc-2
      nset=nset+1
      nset1=nset1+1
      x(nset)=(x(nset1)+x(nset1+1))/2
      y(nset)=(y(nset1)+y(nset1+1))/2
      z(nset)=j*dh
   enddo
   nset=nset+1
   nset1=nset1+1
   x(nset)=(x(nset1)+x(nset1-(nc-1)))/2
   y(nset)=(y(nset1)+y(nset1-(nc-1)))/2
   z(nset)=j*dh
enddo
ntei=nset1 !底辺の中間節点数
!---斜辺の中間節点の座標
nset1=0
do j=0,nh-1
   if (mod(j,2)==0) then
      do i=0,nc-2
         nset=nset+1
         nset1=nset1+1
         x(nset)=(x(nset1)+x(nset1+nc))/2
         y(nset)=(y(nset1)+y(nset1+nc))/2
         z(nset)=j*dh+dh/2
         nset=nset+1
         x(nset)=(x(nset1)+x(nset1+nc+1))/2
         y(nset)=(y(nset1)+y(nset1+nc+1))/2
         z(nset)=j*dh+dh/2
      enddo
      nset=nset+1
      nset1=nset1+1
      x(nset)=(x(nset1)+x(nset1+nc))/2
      y(nset)=(y(nset1)+y(nset1+nc))/2
      z(nset)=j*dh+dh/2
      nset=nset+1
      x(nset)=(x(nset1)+x(nset1+1))/2
      y(nset)=(y(nset1)+y(nset1+1))/2
      z(nset)=j*dh+dh/2
   else
      i=0
      nset=nset+1
      nset1=nset1+1
      x(nset)=(x(nset1)+x(nset1+nc*2-1))/2
      y(nset)=(y(nset1)+y(nset1+nc*2-1))/2
      z(nset)=j*dh+dh/2
      nset=nset+1
      x(nset)=(x(nset1)+x(nset1+nc))/2
      y(nset)=(y(nset1)+y(nset1+nc))/2
      z(nset)=j*dh+dh/2
      do i=1,nc-1
         nset=nset+1
         nset1=nset1+1
         x(nset)=(x(nset1)+x(nset1+nc-1))/2
         y(nset)=(y(nset1)+y(nset1+nc-1))/2
         z(nset)=j*dh+dh/2
         nset=nset+1
         x(nset)=(x(nset1)+x(nset1+nc))/2
         y(nset)=(y(nset1)+y(nset1+nc))/2
         z(nset)=j*dh+dh/2
      enddo
   endif
enddo
nsya=nset1 !斜辺の中間節点数
!---節点番号と座標の出力
print'("*NODE, NSET=Nall")'
do i=1,nset
   print'(i6,a,f17.9,a,f17.9,a,f17.9)',i,',',x(i),',',y(i),',',z(i)
end do
!---要素を構成する節点
nyou=0
do j=0,nh-1
   if (mod(j,2)==0) then
      do i=0,nc-2
         nyou=nyou+1
         mesh(nyou,1)=1+i+nc*j
         mesh(nyou,2)=mesh(nyou,1)+1
         mesh(nyou,3)=mesh(nyou,2)+nc
         mesh(nyou,4)=mesh(nyou,1)+ntyo
         mesh(nyou,5)=mesh(nyou,1)*2+1+ntyo+ntei
         mesh(nyou,6)=mesh(nyou,1)*2+ntyo+ntei
         nyou=nyou+1
         mesh(nyou,1)=2+i+nc*(j+1)
         mesh(nyou,2)=mesh(nyou,1)-1
         mesh(nyou,3)=mesh(nyou,2)-nc
         mesh(nyou,4)=mesh(nyou,3)+ntyo+nc
         mesh(nyou,5)=mesh(nyou,3)*2-1+ntyo+ntei
         mesh(nyou,6)=mesh(nyou,3)*2+ntyo+ntei
      enddo
      i=nc-1
      nyou=nyou+1
      mesh(nyou,1)=1+i+nc*j
      mesh(nyou,2)=mesh(nyou,1)-nc+1
      mesh(nyou,3)=mesh(nyou,1)+1
      mesh(nyou,4)=mesh(nyou,1)+ntyo
      mesh(nyou,5)=mesh(nyou,1)*2+1+ntyo+ntei-2*nc
      mesh(nyou,6)=mesh(nyou,1)*2+ntyo+ntei
      nyou=nyou+1
      mesh(nyou,1)=2+i+nc*j
      mesh(nyou,2)=mesh(nyou,1)+nc-1
      mesh(nyou,3)=mesh(nyou,1)-1
      mesh(nyou,4)=mesh(nyou,3)+ntyo+nc
      mesh(nyou,5)=mesh(nyou,3)*2-1+ntyo+ntei
      mesh(nyou,6)=mesh(nyou,3)*2+ntyo+ntei
   else
      do i=1,nc-1
         nyou=nyou+1
         mesh(nyou,1)=1+i+nc*j
         mesh(nyou,2)=mesh(nyou,1)+nc-1
         mesh(nyou,3)=mesh(nyou,1)-1
         mesh(nyou,4)=mesh(nyou,1)*2+ntyo+ntei-1
         mesh(nyou,5)=mesh(nyou,1)*2+ntyo+ntei-2
         mesh(nyou,6)=mesh(nyou,1)+ntyo-1
         nyou=nyou+1
         mesh(nyou,1)=1+i+nc*(j+1)
         mesh(nyou,2)=mesh(nyou,1)-1
         mesh(nyou,3)=mesh(nyou,1)-nc
         mesh(nyou,4)=mesh(nyou,3)+ntyo+nc-1
         mesh(nyou,5)=mesh(nyou,3)*2+ntyo+ntei-1
         mesh(nyou,6)=mesh(nyou,3)*2+ntyo+ntei
      enddo
      i=0
      nyou=nyou+1
      mesh(nyou,1)=1+i+nc*j
      mesh(nyou,2)=mesh(nyou,1)+nc*2-1
      mesh(nyou,3)=mesh(nyou,1)+nc-1
      mesh(nyou,4)=mesh(nyou,1)*2+ntyo+ntei-1
      mesh(nyou,5)=mesh(nyou,1)*2+ntyo+ntei+(nc-1)*2
      mesh(nyou,6)=mesh(nyou,1)+ntyo+nc-1
      nyou=nyou+1
      mesh(nyou,1)=1+i+nc*(j+1)
      mesh(nyou,2)=mesh(nyou,1)+nc-1
      mesh(nyou,3)=mesh(nyou,1)-nc
      mesh(nyou,4)=mesh(nyou,3)+ntyo+nc*2-1
      mesh(nyou,5)=mesh(nyou,3)*2+ntyo+ntei-1
      mesh(nyou,6)=mesh(nyou,3)*2+ntyo+ntei
   endif
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,nc
   print'(i8,a)',i
enddo
do i=1,nc
   print'(i8,a)',i+ntyo
enddo
!#--載荷点の節点群を出力
print'("*NSET,NSET=saika")'
do i=ntyo-nc+1,ntyo
   print'(i7,a)',i
enddo
do i=ntyo+ntei-nc+1,ntyo+ntei!中間節点にも載荷する場合
   print'(i7,a)',i
enddo
!print'("*RIGID BODY,NSET=saika")'!載荷部分は平らに変形させる(強制変位では使えない)
!---段ごとのz変位を均等にする
print'("*EQUATION")'
do j=1,nh-1
   do i=1,nc
      print'("2")'!項が２つという意味
      print'(i7,A,f10.7,a,i7,A,f10.7)'&
           &,i+nc*j,',3,',real(nh)/real(j),',',ntyo-nc+i,',3,',-1.d0
   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/(nc*2))*i/step!中間節点にかけない(~*1)orかける(~*2)
   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
