!板の非線形座屈(CCX)
! S8

! 実行例：
!  $ ./hitaz<itaz.dat>hitaz.inp

implicit none
!#--入力
real    :: b=100.d-3            !  m :幅
real    :: a=100.d-3            !  m :長さ
real    :: t=1.d-3              !  m :厚さ
integer :: nb=100               !    :幅方向分割数
integer :: na=100               !    :長さ方向分割数
real    :: pcr=0.7385436d-2/0.8 ! MN :座屈荷重(やや大きめに)
real    :: husei=2000           !    :線形座屈の節点変位をこれで割る
real    :: dp=1.d-3             !    :1を全体とした時の荷重増分?
integer :: nstep=100            !    :出力ステップ
!------
real    :: x(1000000),y(1000000),z(1000000)
real    :: u(1000000),v(1000000),w(1000000)
integer :: mesh(1000000,8)
real    :: db,da,ts
integer :: i,j,k,ndan,nset,nyou,gomi

db=b/nb !  m :一要素の幅
da=a/na !  m :一要素の長さ

!#--datファイルの読み込み
do i=1,30421
   read(*,'()')
enddo
do i=1,(nb*3+2)*na+(nb*2+1)
   read(*,*) gomi,u(i),v(i),w(i)
enddo

!---節点座標の計算
nset=0
do j=0,na*2
   if (mod(j,2)==1) then
      do i=0,nb
         nset=nset+1
         x(nset)=db*i+u(nset)/husei
         y(nset)=v(nset)/husei
         z(nset)=da/2*j+w(nset)/husei
      enddo
   else
      do i=0,nb*2
         nset=nset+1
         x(nset)=db/2*i+u(nset)/husei
         y(nset)=v(nset)/husei
         z(nset)=da/2*j+w(nset)/husei
      enddo
   endif
enddo
!---節点の出力
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)
enddo
!---要素を構成する節点
ndan=nb*3+2
nyou=0
do j=0,na-1
   do i=0,nb-1
      nyou=nyou+1
      mesh(nyou,1)=1+i*2+ndan*j
      mesh(nyou,2)=mesh(nyou,1)+2
      mesh(nyou,3)=mesh(nyou,2)+ndan
      mesh(nyou,4)=mesh(nyou,3)-2
      mesh(nyou,5)=mesh(nyou,1)+1
      mesh(nyou,6)=2+i+(nb*2+1)+ndan*j
      mesh(nyou,7)=mesh(nyou,3)-1
      mesh(nyou,8)=mesh(nyou,6)-1
   enddo
enddo
!---要素の出力
print'("*ELEMENT, TYPE=S8, 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),',',&
       &mesh(i,7),',',&
       &mesh(i,8)
end do
!---境界条件の節点群を出力
print'("*NSET,NSET=kotei1")'
do i=1,nb*2+1
   print'(i8,a)',i
enddo
print'("*NSET,NSET=kotei2")'
do i=0,na
   print'(i8,a)',ndan*i+1
enddo
do i=0,na-1
   print'(i8,a)',ndan*i+(nb*2+1)+1
enddo
print'("*NSET,NSET=kotei3")'
do i=0,na
   print'(i8,a)',ndan*i+nb*2+1
enddo
do i=0,na-1
   print'(i8,a)',ndan*i+ndan
enddo
!---荷重/境界条件の節点群を出力
print'("*NSET,NSET=saika")'
do i=nset-nb*2,nset
   print'(i8,a)',i
enddo
!---境界条件
print'("*BOUNDARY")'
print'("kotei1,2,3")'
print'("kotei2,1,2")'
print'("kotei3,2")'
print'("saika,2")'
print'("*EQUATION")'!kotei3のx方向変位、saikaのz方向変位を等しく
do i=nb*2+1,ndan*(na-1)+(nb*2+1),ndan!kotei3
   print'("2")'
   print'(i6,A,i6,A)',i,',1,1,',i+nb+1,',1,-1'
   print'("2")'
   print'(i6,A,i6,A)',i+nb+1,',1,1,',i+nb*3+2,',1,-1'
enddo
do i=nset-nb*2,nset-1!saika
   print'("2")'
   print'(i6,A,i6,A)',i,',3,1,',i+1,',3,-1'
enddo
!#---材料定数と厚さ
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206d3, 0.3d0")'
print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL1")'
print'(1pe13.5)',t
!---出力ステップ
print'("*TIME POINTS,NAME=T1")'
do i=1,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo
!---載荷条件
print'("*STEP,NLGEOM")'
print'("*STATIC,SOLVER=SPOOLES")'
!print'(1pe13.5,a,1pe13.5)',dp,',',1.0
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',-pcr/(nb*2+1)
!---出力の設定
print'("*NODE PRINT,NSET=Nall,TIMEPOINTS=T1")'
print'("U")'
!print'("*EL PRINT,ELSET=Eall,TIMEPOINTS=T1")'
!print'("S")'
print'("*NODE FILE,TIMEPOINTS=T1")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T1")'
print'("S")'
print'("*END STEP")'
end
