!CLTの静的載荷の1/4解析(CCX)の幅を変化させて、一気に計算。
! 載荷幅と支間距離は実験に合わせたが、固定面は途中で途切れない。
! 中央集中荷重だったものを分布荷重に変えた。

! 実行例：
!  $ gf clthaba clthaba.f90
!  $ ./clthaba
!  $ chmod 744 job
!  $ ./job
implicit none
!#--入力
real    :: b1s=200.d-3/2    ! m      :短い幅/2(x)
real    :: b1l=4000.d-3/2   ! m      :長い幅/2(x)
integer :: nb1=19           !        :長さ変化の分割数
real    :: h=240.d-3        ! m      :全層の高さ(y)
real    :: ell=2010.d-3/2   ! m      :短い長さ/2(z)
real    :: b2=3380.d-3      ! m      :固定しているH鋼の長さ(x)
real    :: ellb=80.d-3      ! m      :端からH鋼に固定しているボルトまでの距離(z)
real    :: ellh=300.d-3/2   ! m      :H鋼のフランジ幅/2(z)
integer :: nsou=8           !        :層数
integer :: nx=100           !        :幅方向分割数
integer :: ny=3             !        :1層ごとの高さ方向分割数
!integer :: ny=24            !        :1層ごとの高さ方向分割数
integer :: nz=100           !        :長さ方向分割数
real    :: q=1.0d-3         ! MN/m^2 :荷重(活荷重のイメージで)
real    :: hennix=200.d-3/2 ! m      :変位を見たい箇所の幅/2(x)
real    :: henniz=500.d-3/2 ! m      :変位を見たい箇所の幅/2(z)
real    :: exx=5.0d3/25     ! MPa    :x軸(幅)方向のヤング率
!real    :: exx=1.7d3        ! MPa    :x軸(幅)方向のヤング率
real    :: eyy=5.0d3/25     ! MPa    :y軸(高さ)方向のヤング率
real    :: ezz=5.0d3        ! MPa    :z軸(長さ)方向のヤング率
!real    :: ezz=3.5d3        ! MPa    :z軸(長さ)方向のヤング率
real    :: poizx=0.4d0      !        :ポアソン比
!---
real poixy,poixz,poiyx,poiyz,poizy,gxy,gxz,gyz
real :: x(1000000),y(1000000),z(1000000) !節点座標
integer :: mesh(1000000,8)               !要素と要素を構成する節点
real :: e66(6,6),b(3,3),dx,dy,dz,b1,dell
integer i,j,k,l,m,nset,nyou,nsaika
character :: filename1*10,filename2*30
!---ポアソン比
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poizx*eyy/ezz
poiyz=poiyx
!---せん断弾性係数
gxy=eyy*25.d0/15
gxz=gxy
gyz=gxy
!---出力ファイルを初期化
open(8,file='job',status='replace')
close(8)
!---幅を変化させるループ
do m=0,nb1
   b1=b1s+(b1l-b1s)/nb1*m
   !#--出力ファイル名
   write(filename1,'(i10)') int(b1*1.d3)
   open(7,file="haba"//trim(adjustl(filename1))//".inp", status='replace')
   !---要素長さ
   dx=b1/nx
   dy=h/(ny*nsou)
   dz=ell/nz
   !---節点座標の計算
   nset=0
   do k=0,nz
      do j=0,ny*nsou
         do i=0,nx
            nset=1+nset
            x(nset)=i*dx
            y(nset)=j*dy
            z(nset)=k*dz
         enddo
      enddo
   enddo
   !---要素を構成する節点の計算
   nyou=0
   do k=0,nz-1
      do j=0,ny*nsou-1
         do i=1,nx
            nyou=1+nyou
            mesh(nyou,1)=i+j*(nx+1)+k*((nx+1)*(ny*nsou+1))
            mesh(nyou,2)=mesh(nyou,1)+(nx+1)
            mesh(nyou,3)=mesh(nyou,2)+(nx+1)*(ny*nsou+1)
            mesh(nyou,4)=mesh(nyou,1)+(nx+1)*(ny*nsou+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*nsou+1)
            mesh(nyou,8)=mesh(nyou,5)+(nx+1)*(ny*nsou+1)
         enddo
      enddo
   enddo
   !---節点座標の出力
   write(7,'("*NODE, NSET=Nall")')
   do i=1,nset
      write(7,'(i6,a,f9.6,a,f9.6,a,f9.6)'),&
          &i,',',x(i),',',y(i),',',z(i)
   end do
   !---要素の出力
   do k=0,nsou-1
      write(7,'(A,i1)'),"*ELEMENT, TYPE=C3D8, ELSET=E",k+1
      do j=0,nz-1
         do i=1,nx*ny
            l=i+j*(nx*ny*nsou)+k*(nx*ny)
            write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)'),&
                &l,',',&
                &mesh(l,1),',',&
                &mesh(l,2),',',&
                &mesh(l,3),',',&
                &mesh(l,4),',',&
                &mesh(l,5),',',&
                &mesh(l,6),',',&
                &mesh(l,7),',',&
                &mesh(l,8)
         enddo
      enddo
   enddo
   !---境界条件の節点群
   write(7,'("*NSET,NSET=kotei1")')!H鋼の半面
   do i=1,nset
      if (y(i)>h-dy/3.and.y(i)<h+dy/3.and.&
          &z(i)>ell-ellh-dz/3.and.z(i)<ell-ellb+dz/3) then
         write(7,'(i6,a)'),i
      endif
   end do
   write(7,'("*NSET,NSET=kotei2")')!x対称面
   do i=1,nset
      if (x(i)>-dx/3.and.x(i)<dx/3) then
         write(7,'(i6,a)'),i
      endif
   end do
   write(7,'("*NSET,NSET=kotei3")')!z対称面
   do i=1,nset
      if (z(i)>-dz/3.and.z(i)<dz/3) then
         write(7,'(i6,a)'),i
      endif
   end do
   !---載荷条件の節点群
   write(7,'("*NSET,NSET=saika")')
   nsaika=0
   do i=1,nset
      if (y(i)>-dy/3.and.y(i)<dy/3) then
         nsaika=nsaika+1
         write(7,'(i6,a)'),i
      endif
   end do
   !---変位を見たい箇所の節点群
   write(7,'("*NSET,NSET=henni")')
   do i=1,nset
      if (y(i)>-dy/3.and.y(i)<dy/3.and.&
          &z(i)>=-dz/3.and.z(i)<=henniz+dz/3.and.&
          &x(i)>=-dx/3.and.x(i)<=hennix+dx/3) then
         write(7,'(i6,a)'),i
      endif
   end do
   !---境界条件
   write(7,'("*BOUNDARY")')
   write(7,'("kotei1,1,2")')
   write(7,'("kotei2,1")')
   write(7,'("kotei3,3")')
   !---材料設定
   write(7,'("*MATERIAL,NAME=EL1")')
   !write(7,'("*ELASTIC")')!等方性
   !write(7,'(f12.4,a,f12.4)'),ezz,',',poizx
   write(7,'("*ELASTIC,TYPE=ORTHO")')
   call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
   & gxz,gxy,poiyx,poizx,poizy)
   write(7,'(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),','
   write(7,'(f12.4)'),e66(4,4)
   !#--局所座標系の設定
   write(7,'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")')
   write(7,'("1.,0.,0.,0.,1.,0.")')!全体座標系と一致
   write(7,'("*ORIENTATION,NAME=OR2,SYSTEM=RECTANGULAR")')
   write(7,'("0.,0.,-1.,0.,1.,0.")')!y軸回りに90度回転
   write(7,'("*ORIENTATION,NAME=OR3,SYSTEM=RECTANGULAR")')
   write(7,'("1.,0.,0.,0.,0.,-1.")')!x軸回りに90度回転
   write(7,'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1,ORIENTATION=OR1")')
   write(7,'("*SOLID SECTION,ELSET=E2,MATERIAL=EL1,ORIENTATION=OR2")')
   write(7,'("*SOLID SECTION,ELSET=E3,MATERIAL=EL1,ORIENTATION=OR1")')
   write(7,'("*SOLID SECTION,ELSET=E4,MATERIAL=EL1,ORIENTATION=OR2")')
   write(7,'("*SOLID SECTION,ELSET=E5,MATERIAL=EL1,ORIENTATION=OR2")')
   write(7,'("*SOLID SECTION,ELSET=E6,MATERIAL=EL1,ORIENTATION=OR1")')
   write(7,'("*SOLID SECTION,ELSET=E7,MATERIAL=EL1,ORIENTATION=OR2")')
   write(7,'("*SOLID SECTION,ELSET=E8,MATERIAL=EL1,ORIENTATION=OR1")')
   !---載荷条件
   write(7,'("*STEP")')
   write(7,'("*STATIC,SOLVER=SPOOLES")')
   write(7,'("*CLOAD")')
   write(7,'(A,a,i6,a,1pe13.5)'),'saika',',',2,',',q*ell*b1/nsaika
   !#--datファイルの設定
   write(7,'("*NODE PRINT,NSET=henni")')
   write(7,'("U")')
   !write(7,'("*EL PRINT,ELSET=")')
   !write(7,'("S")')
   write(7,'("*NODE FILE")')
   write(7,'("U")')
   write(7,'("*EL FILE")')
   write(7,'("S")')
   write(7,'("*END STEP")')
   !#--jobの名前
   open(8,file='job',access='append')
   filename2="ccx haba"//trim(adjustl(filename1))&
       &//">haba"//trim(adjustl(filename1))//".kekka"
   write(8,'(A)'),filename2
   close(8)
   close(7)
enddo!幅変化のdoループ
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
