!座屈の非線形解析(CCX)
! 線形座屈解析での hoge.inp と hoge.dat から、
! 初期たわみを与えたhoge2.inpをつくるプログラム。
! 多分、要素は何でもいける。

!使い方:
! 0. 線形座屈解析をするとき、全節点変位をdatファイルに出力するようにする。
! 1. $ vi syokitawami.f90
!    hoge.inpとhoge.datをよく見ながら、syokihusei.f90の#のところを書き換える。
! 2. $ gfortran syokitawami.f90 -o syokitawami
! 3. $ ./syokitawami>hoge2.inp
! 4.
! 4.1 $ cgx -c hoge2.inp
! 4.2 "plot n yomi"と打ち、変位を読み取りたい座屈箇所の節点が表示されるかを確認しておく。
!     節点が選択されていなかったり、あるいは多すぎたりした場合は、1. に戻る。
! 5. $ vi hoge2.inp
! 5.1. "*STEP"を"*STEP,NLGEOM,INC=1000"(INCはnstepより多ければいい)に書き換えて、
!      必要に応じてその下の行に荷重増分を追加する。
! 5.2. "*BUCKLE"の2行を消す。
! 5.3. 元の"*CLOAD"を"*END STEP"以下に書き込んだものに置き換える。
! 5.4. "*END STEP"以下に書き込んだ出力ステップ数の設定を"*STEP"の上に切り取って貼り付ける。
! 5.5. "*NODE PRINT,NSET=saika"を"NODE PRINT,NSET=yomi,TIMEPOINTS=T1"に書き換える。
! 6. $ ccx hoge2

implicit real*8(a-h,o-z)
dimension yomi(100)
character(200) itigyou
!#--手入力
nnode=585843         !   :節点数
nsaika=57513         !   :載荷節点数
husei=2000           !   :線形座屈の節点変位をこれで割る
ngyou=1171712        !   :datファイルの始めのいらない行数(モードによって変える)
!zakutu=0.2305162d0   ! MN:見たいモードの座屈荷重
!zakutu=0.2305162d0/4 ! MN:見たいモードの座屈荷重(1/4解析)
yuki=0.1750624d3     !  m:積雪深
a=0.42d0*3.65d0      !m^2:上面の面積
p=zakutu/0.75        ! MN:集中荷重.座屈後挙動も見たいのでちょっと大きめに荷重を与える
q=3.5d-3*a*yuki/0.75 !MN/m^2:等分布荷重.座屈後挙動も見たいのでちょっと大きめに荷重を与える
nstep=200            !   :出力ステップ数
!#--読み込みファイル名
open(11,file='kou30_yuki.inp',status='old')!読み取りのみ可
open(12,file='kou30_yuki.dat',status='old')
!---datファイルのデータの読み飛ばし
do i=1,ngyou
   read(12,'()')
enddo
!---元の座標と節点変位の読み込み
read(11,'()')
print'("*NODE, NSET=Nall")'
nyomi=0
do i=1,nnode
   read(12,*) gomi,u,v,w!節点変位
   read(11,*) gomi,x,y,z!元の座標
   !#--datファイルの読み取り点をつくる
   if (x>-1.d-3.and.x<+1.d-3.and.&
      &y>0.25-1.d-3.and.y<0.25+1.d-3.and.&
      &z>0.73-1.d-3.and.z<0.73+1.d-3) then
      nyomi=nyomi+1
      yomi(nyomi)=i
   endif
   !#--節点座標の出力
   print'(i6,a,f17.9,a,f17.9,a,f17.9)',&
   !    &i,',',x+u/husei,',',y+v/husei,',',z+w/husei
       &i,',',x-u/husei,',',y-v/husei,',',z-w/husei!反対方向に初期たわみを与える場合
enddo
!---読み取り点の出力
print'("*NSET, NSET=yomi")'
do i=1,nyomi
   print'(i6,a)',int(yomi(i)),','
enddo
!---元inpファイルの節点情報以外の全てを新inpファイルに書き写し
do!無限ループ
   read(11,'(a)',iostat=ios) itigyou!一行そっくりそのまま読み取る
   if (ios<0) exit!無限ループから抜ける
      print*,trim(adjustl(itigyou))!空白を取り除いてから、書き込む
enddo
!#--荷重の設定
print'("*CLOAD")'
!print'(A,a,i6,a,1pe17.9)','sensaika',',',2,',',p/nsaika
print'(A,a,i6,a,1pe17.9)','mensaika',',',2,',',q/nsaika
!---出力ステップの設定
print'("*TIME POINTS,NAME=T1")'
do i=0,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo

close(11)
close(12)
end

