!片持ち梁の横ねじれ座屈の非線形解析(CCX)
! 弾性体、等方性材料

! 実行例：
! yoko.datのデータを入力しておいて、コンパイル
! $ ./hyoko<yoko.dat>hyoko.inp

implicit real*8(a-h,o-z)
dimension x(1000000),y(1000000),z(1000000),e66(6,6),b(3,3)&
&,mesh(30000000,8),u(1000000),v(1000000),w(1000000)
!プログラムはメートル、MPa、
zakutu=0.4682880E-01!MN:起こしたいモードの座屈荷重
!zakutu=0.1212088E+00!MN:起こしたいモードの座屈荷重
!zakutu=0.2000032E+00!MN:起こしたいモードの座屈荷重
!zakutu=0.2845026E+00!MN:起こしたいモードの座屈荷重
!zakutu=0.3756069E+00!MN:起こしたいモードの座屈荷重
zoubun=0.005        !1を全体とした時の荷重増分?
p=zakutu/0.8        !MN:ちょっと大きめに荷重をかける
husei=5000          !線形座屈の節点変位をこれで割る
nstep=100            !出力ステップ数
!寸法
b1=2.0d-2           !幅(m)
h=25.d-2            !高さ(m)
ell=2.0d0           !長さ(m)
!xyz方向の分割数
nx=8
ny=30
nz=200
!datファイルのデータの読み飛ばし
do i=1,56099 !初めの読み飛ばす行数(モードによって変わる)
read(*,'()')
enddo
!datファイルの節点変位の読み込み
do i=1,(nx+1)*(ny+1)*(nz+1) !節点数
read(*,*) a,u(i),v(i),w(i)
enddo
!1つあたりの要素の長さ
xn=b1/nx
yn=h/ny
zn=ell/nz
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
do k=0,nz!z
do j=0,ny!
do i=0,nx!x方向分割数を担当
nset=1+i+j*(nx+1)+k*((nx+1)*(ny+1))
x(nset)=i*xn+u(nset)/husei
y(nset)=j*yn+v(nset)/husei
z(nset)=k*zn+w(nset)/husei
end do
end do
end do 
!i5~6が整数,aが文字,f9,6が小数.1peとかはEの＊乗
!接点番号と座標の出力
nsettenzenbu=(nx+1)*(ny+1)*(nz+1)
print'("*NODE, NSET=Nall")'
do i=1,nsettenzenbu
print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
    &i,',',x(i),',',y(i),',',z(i)
nsettensuu=i
end do
!要素番号とその要素の8接点の節点番号
nyou=0
do k=0,nz-1
do j=0,ny-1
do i=1,nx
nyou=nyou+1
mesh(nyou,1)=i+j*(nx+1)+k*((nx+1)*(ny+1))
mesh(nyou,2)=mesh(nyou,1)+(nx+1)
mesh(nyou,3)=mesh(nyou,2)+(nx+1)*(ny+1)
mesh(nyou,4)=mesh(nyou,1)+(nx+1)*(ny+1)
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nx+1
mesh(nyou,7)=mesh(nyou,6)+(nx+1)*(ny+1)
mesh(nyou,8)=mesh(nyou,5)+(nx+1)*(ny+1)
end do
end do
end do
!要素の出力
print'("*ELEMENT, TYPE=C3D8, ELSET=E1")'
do i=1,nyou
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
nitimen=(nx+1)*(ny+1)
!固定端の節点群を書きだす
print'("*NSET,NSET=kotei")'
do i=1,nitimen
nkotei=i
print'(i6,a)',nkotei
end do
!載荷点の節点群を書きだす
print'("*NSET,NSET=saika")'
do i=1,nitimen
nsaika=nitimen*nz+i
print'(i6,a)',nsaika
end do
!境界条件
print'("*BOUNDARY")'
print'("kotei,1,3")'
! 木材の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")'
!出力ステップ
print'("*TIME POINTS,NAME=T1")'
do i=0,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo
! 載荷条件
print'("*STEP,NLGEOM")'
print'("*STATIC,SOLVER=SPOOLES")'
!print'(1pe13.5,a,1pe13.5)',zoubun,',',1.0
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nitimen
! datのオプション
print'("*NODE PRINT,NSET=saika,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,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'
!
end

