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

! 実行例：
! szmoku.datのデータを入力しておいて、コンパイル
! $ ./hzmoku<szmoku.dat>hzmoku.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.7586853d-3 !MN:起こしたいモードの座屈荷重
!zakutu=0.2038648d-2 !MN:起こしたいモードの座屈荷重
!zakutu=0.3523796d-2 !MN:起こしたいモードの座屈荷重
!zakutu=0.5226363d-2 !MN:起こしたいモードの座屈荷重
!zakutu=0.7132325d-2 !MN:起こしたいモードの座屈荷重
zoubun=0.01         !1を全体とした時の荷重増分?
p=zakutu/0.8        !MN:ちょっと大きめに荷重をかける
husei=5000          !線形座屈の節点変位をこれで割る
!寸法
b1=2.0d-2           !幅(m)
h=25.d-2            !高さ(m)
ell=2.0d0           !長さ(m)
!ヤング率
ezz=7.694d3
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!xyz方向の要素分割数
nx=6
ny=30
nz=200
!datファイルのデータの読み飛ばし
do i=1,43637 !初めの読み飛ばす行数(モードによって変わる)
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,TYPE=ORTHO")'
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
print'(f12.4,a,f12.4,a,f12.4,a,&
&f12.4,a,f12.4,a,f12.4,a,f12.4,&
& a ,f12.4,a)',&
&e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
&,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
print'(f12.4)',e66(4,4)
print'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")'
! 載荷条件
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")'
print'("U")'
print'("*EL PRINT,ELSET=Eall")'
print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'
!
end

!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz

!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
&+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
&-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
return
end
