!弾性床上の梁の座屈(CCX)
! 床:木材 / 梁:鋼材 を想定
implicit none
!#--入力
real    :: b1=200.d-4      ! m   :床の幅(x)
real    :: ell1=400.d-3    ! m   :床の長さ(z)
real    :: h1=120.d-4      ! m   :床の高さ(y)
real    :: b2=120.d-4      ! m   :梁の幅(x)
real    :: ell2=315.d-3    ! m   :梁の長さ(z)
real    :: h2=9.d-3        ! m   :梁の高さ(y)
integer :: nx1=5           !     :床の幅方向分割数
integer :: ny1=4           !     :床の高さ方向分割数
integer :: nz1=100         !     :床の長さ方向分割数
integer :: nx2=5           !     :梁の幅方向分割数
integer :: ny2=4           ! 偶数:梁の高さ方向分割数
integer :: nz2=100         !     :梁の長さ方向分割数
real    :: p=20.d-3          ! MN  :荷重
integer :: nstep=200      !     :出力ステップ数
real    :: syokix=100.d-4  ! m   :梁の初期中心位置x
real    :: syokiz=200.d-3  ! m   :梁の初期中心位置z
real    :: exx=6.0d3/25    ! MPa :x軸方向のヤング率
real    :: eyy=6.0d3/25    ! MPa :y軸方向のヤング率
real    :: ezz=6.0d3       ! MPa :z軸方向のヤング率
real    :: poizx=0.4d0     !     :ポアソン比
real    :: ekou=205.d3     ! MPa :鋼板のヤング率
real    :: poikou=0.3d0    !     :鋼板のポアソン比
!---
real poixy,poixz,poiyx,poiyz,poizy,gxy,gxz,gyz
real :: x1(1000000),y1(1000000),z1(1000000) !相手材の節点座標
real :: x2(1000000),y2(1000000),z2(1000000) !相手材の節点座標
integer :: mesh1(1000000,8),mesh2(1000000,8) !要素と要素を構成する節点
real :: e66(6,6),b(3,3),dx1,dy1,dz1,dx2,dy2,dz2,ts
integer i,j,k,l,nset1,nset2,nyou1,nyou2,nsaika
!---ポアソン比
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poixy
poiyz=poixy
!---せん断弾性係数
gxy=ezz/15
gxz=gxy
gyz=gxy
!---要素長さ
dx1=b1/nx1
dy1=h1/ny1
dz1=ell1/nz1
dx2=b2/nx2
dy2=h2/ny2
dz2=ell2/nz2
!---節点座標の計算
nset1=0
do k=0,nz1
   do j=0,ny1
      do i=0,nx1
         nset1=1+nset1
         x1(nset1)=i*dx1
         y1(nset1)=j*dy1
         z1(nset1)=k*dz1
      enddo
   enddo
enddo
nset2=0
do k=0,nz2
   do j=0,ny2
      do i=0,nx2
         nset2=1+nset2
         x2(nset2)=syokix-b2/2+i*dx2
         y2(nset2)=h1+j*dy2
         z2(nset2)=syokiz-ell2/2+k*dz2
      enddo
   enddo
enddo
!---要素を構成する節点の計算
nyou1=0
do k=0,nz1-1
   do j=0,ny1-1
      do i=1,nx1
         nyou1=1+nyou1
         mesh1(nyou1,1)=i+j*(nx1+1)+k*((nx1+1)*(ny1+1))
         mesh1(nyou1,2)=mesh1(nyou1,1)+(nx1+1)
         mesh1(nyou1,3)=mesh1(nyou1,2)+(nx1+1)*(ny1+1)
         mesh1(nyou1,4)=mesh1(nyou1,1)+(nx1+1)*(ny1+1)
         mesh1(nyou1,5)=mesh1(nyou1,1)+1
         mesh1(nyou1,6)=mesh1(nyou1,5)+nx1+1
         mesh1(nyou1,7)=mesh1(nyou1,6)+(nx1+1)*(ny1+1)
         mesh1(nyou1,8)=mesh1(nyou1,5)+(nx1+1)*(ny1+1)
      enddo
   enddo
enddo
nyou2=0
do k=0,nz2-1
   do j=0,ny2-1
      do i=1,nx2
         nyou2=1+nyou2
         mesh2(nyou2,1)=nset1+i+j*(nx2+1)+k*((nx2+1)*(ny2+1))
         mesh2(nyou2,2)=mesh2(nyou2,1)+(nx2+1)
         mesh2(nyou2,3)=mesh2(nyou2,2)+(nx2+1)*(ny2+1)
         mesh2(nyou2,4)=mesh2(nyou2,1)+(nx2+1)*(ny2+1)
         mesh2(nyou2,5)=mesh2(nyou2,1)+1
         mesh2(nyou2,6)=mesh2(nyou2,5)+nx2+1
         mesh2(nyou2,7)=mesh2(nyou2,6)+(nx2+1)*(ny2+1)
         mesh2(nyou2,8)=mesh2(nyou2,5)+(nx2+1)*(ny2+1)
      enddo
   enddo
enddo
!---節点座標の出力
print'("*NODE, NSET=N1")'
do i=1,nset1
   print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
       &i,',',x1(i),',',y1(i),',',z1(i)
end do
print'("*NODE, NSET=N2")'
do i=1,nset2
   print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
       &nset1+i,',',x2(i),',',y2(i),',',z2(i)
end do
!---要素の出力
print'(A)',"*ELEMENT, TYPE=C3D8, ELSET=E1"
do i=1,nyou1
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
      &i,',',&
      &mesh1(i,1),',',&
      &mesh1(i,2),',',&
      &mesh1(i,3),',',&
      &mesh1(i,4),',',&
      &mesh1(i,5),',',&
      &mesh1(i,6),',',&
      &mesh1(i,7),',',&
      &mesh1(i,8)
enddo
print'(A)',"*ELEMENT, TYPE=C3D8, ELSET=E2"
do i=1,nyou2
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
      &nyou1+i,',',&
      &mesh2(i,1),',',&
      &mesh2(i,2),',',&
      &mesh2(i,3),',',&
      &mesh2(i,4),',',&
      &mesh2(i,5),',',&
      &mesh2(i,6),',',&
      &mesh2(i,7),',',&
      &mesh2(i,8)
enddo
!---マスター面群
print'(A)',"*SURFACE,NAME=SS"
do k=0,nz1-1
   do i=1,nx1
     l=nx1*(ny1-1)+nx1*ny1*k+i
      print'(i6,A)',l,',S4'
   enddo
enddo
!---スレーブ面群
print'(A)',"*SURFACE,NAME=SM"
do k=0,nz2-1
   do i=1,nx2
      l=nyou1+nx2*ny2*k+i
      print'(i6,A)',l,',S6'
   enddo
enddo
!---境界条件の節点群
print'("*NSET,NSET=kotei0")'!床の下
do i=1,nset1
   if (y1(i)>-dy1/3.and.y1(i)<dy1/3) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=kotei1")'!梁の片端
!do i=nset1+1,nset1+(nx2+1)*(ny2+1)!全面を境界条件に
!do i=nset1+1,nset1+nx2+1!床と接している辺を境界条件に
do i=nset1+(nx2+1)*ny2/2+1,nset1+(nx2+1)*ny2/2+nx2+1!中立軸を境界条件に
   print'(i6,a)',i
end do
print'("*NSET,NSET=kotei2")'!梁のもう片端
!do i=nset1+nset2-(nx2+1)*(ny2+1)+1,nset1+nset2-(nx2+1)*ny2!床と接している辺を境界条件に
do i=nset1+nset2-(nx2+1)*(ny2/2+1)+1,nset1+nset2-(nx2+1)*ny2/2!中立軸を境界条件に
   print'(i6,a)',i
end do
!---載荷条件の節点群
print'("*NSET,NSET=saika")'
nsaika=0
do i=1,nset2
   if (z2(i)>syokiz+ell2/2-dz2/3.and.y2(i)<syokiz+ell2/2+dz2/3) then
      nsaika=nsaika+1
      print'(i6,a)',nset1+i
   endif
end do
!---境界条件
print'("*BOUNDARY")'
print'("kotei0,1,3")'!床の下
print'("kotei1,1,3")'!梁の片端
print'("kotei2,1,2")'!梁のもう片端
!#--接触設定
print'("*CONTACT PAIR,INTERACTION=SI1")'
print'("SS,SM")'!スレーブ面,マスター面
print'("*SURFACE INTERACTION,NAME=SI1")'
print'("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=EXPONENTIAL")'
print'("1.e-4,1.e0")'!接触と判定する距離(m),距離0の時の応力?
print'("*FRICTION")'
print'("0.35,1000")'!摩擦係数,謎パラメータ
!---材料設定
print'("*MATERIAL,NAME=ELkou")'!鋼材
print'("*ELASTIC")'
print'(f12.4,a,f12.4)',ekou,',',poikou
print'("*SOLID SECTION,ELSET=E2,MATERIAL=ELkou")'
print'("*MATERIAL,NAME=ELmoku")'!木材
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'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")'
print'("1.,0.,0.,0.,1.,0.")'!全体座標系と一致
print'("*ORIENTATION,NAME=OR2,SYSTEM=RECTANGULAR")'
print'("0.,0.,-1.,0.,1.,0.")'!y軸回りに90度回転
print'("*ORIENTATION,NAME=OR3,SYSTEM=RECTANGULAR")'
print'("1.,0.,0.,0.,0.,-1.")'!x軸回りに90度回転
print'("*SOLID SECTION,ELSET=E1,MATERIAL=ELmoku,ORIENTATION=OR1")'
!---出力ステップ
print'("*TIME POINTS,NAME=T1")'
do i=0,nstep
   ts=real(i)/nstep
   print'(f12.4,a)',ts,','
enddo
!#--載荷条件
print'("*STEP,NLGEOM,INC=10000")'
!print'("*STEP")'
print'("*STATIC")'
!print'("0.01,1.0")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',-p/nsaika
!print'("*BUCKLE")'
!print'("5,0.01")'!物体の大きさによっては、右の値(Accuracy desired)を*10^?に変化させる。
!#--datファイルの設定
!print'("*NODE PRINT,NSET=N2")'
print'("*NODE PRINT,NSET=N2,TIMEPOINTS=T1")'
print'("U")'
!print'("*EL PRINT,ELSET=")'
!print'("*EL PRINT,ELSET=,TIMEPOINTS=T1")'
!print'("S")'
!print'("*NODE FILE")'
print'("*NODE FILE,TIMEPOINTS=T1")'
print'("U")'
!print'("*EL FILE")'
print'("*EL FILE,TIMEPOINTS=T1")'
print'("S")'
print'("*END STEP")'

end

!-------- 直交異方性の定数を出すプログラム --------
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit none
real :: e66(6,6),b(3,3)
real exx,eyy,ezz,poixy,poixz,poiyz,gyz,gxz,gxy,poiyx,poizx,poizy
real det
integer i,j
!implicit real*8(a-h,o-z)
! 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=0.d0
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
