!CLTの海老名1層検算(CCX)
! (若干寸法違うけど)接続金物(板要素)有り、ゴム有り
! 弾塑性への移行等を想定して、半解析に変えた。
! 境界条件をぴったりに合わせるように変えた。
! 変位計の位置で変位を出力するように変えた。(載荷位置は手動で変えるべし)
implicit none
!#--入力
real    :: b1=3380.d-3    ! m   :x軸(橋軸)方向の長さ1(H形鋼の長さ)
real    :: b2=310.d-3     ! m   :x軸(橋軸)方向の長さ2(両あまり)
real    :: h1=240.d-3     ! m   :全層の高さ(y軸)
real    :: h2=20.d-3      ! m   :ゴムの高さ(y軸)
real    :: ell1=860.d-3   ! m   :z軸(支間)方向の長さ1
real    :: ell2=60.d-3    ! m   :z軸(支間)方向の長さ2(ボルト~ゴム)
real    :: ell3=80.d-3    ! m   :z軸(支間)方向の長さ3(端~ボルト)
real    :: t=18.d-3       ! m   :接続金物の厚さ
real    :: b3=120.d-3     ! m   :接続金物の上部の幅
real    :: b4=320.d-3     ! m   :接続金物の下部の幅
integer :: nsou=8         !     :層数
integer :: nx1=50/2       !     :x軸(橋軸)方向分割数1/2
integer :: nx2=3          !     :x軸(橋軸)方向分割数2
integer :: nx4=8          ! 偶数:接続金物上部のx軸(橋軸)方向分割数
integer :: ny1=2          ! 偶数:y軸(高さ)方向分割数
integer :: ny2=2          !     :ゴムのy軸(高さ)方向分割数
integer :: nz1=50         ! 偶数:z軸(支間)方向分割数1
integer :: nz2=2          ! 偶数:z軸(支間)方向分割数2
integer :: nz3=2          ! 偶数:z軸(支間)方向分割数3
real    :: p=8.2d-3*9.8/2 ! MN  :荷重/2
real    :: saikax=200.d-3 ! m   :載荷板の幅(x)
real    :: saikaz=500.d-3 ! m   :載荷板の幅(z)
real    :: itix=2000.d-3  ! m   :橋軸方向(x)の載荷板中心位置
!real    :: itix=1000.d-3  ! m   :橋軸方向(x)の載荷板中心位置
real    :: itiz=0.d-3     ! m   :支間方向(z)の載荷板中心位置
real    :: exx=1.7d3      ! MPa :x軸(橋軸)方向のヤング率
real    :: eyy=5.0d3/25   ! MPa :y軸(高さ)方向のヤング率
real    :: ezz=3.5d3      ! MPa :z軸(支間)方向のヤング率
real    :: poizx=0.4d0    !     :ポアソン比
real    :: ekou=206.d3    ! MPa :鋼材のヤング率
real    :: poikou=0.3d0   !     :鋼材のポアソン比
real    :: egom=1.d0      ! MPa :ゴムのヤング率
real    :: poigom=0.499d0 !     :ゴムのポアソン比
!---
real poixy,poixz,poiyx,poiyz,poizy,gxy,gxz,gyz
real :: x(1000000),y(1000000),z(1000000) !節点座標
integer :: mesh(1000000,8)               !CLTの要素と要素を構成する節点
integer :: meshk(10000,6)                !CLT角の要素と要素を構成する節点
integer :: meshj(10000,8)                !接続金具の要素と要素を構成する節点
integer :: meshg(10000,8)                !ゴムの要素と要素を構成する節点
real :: e66(6,6),b(3,3),dx1,dx2,dx3,dx4,dx5,dy1,dy2,dz1,dz2,dz3,kaku
integer i,j,k,l,j1,i2,nx,nx3,nz
integer nset,nsetclt,nsetgom,nyou,nyouclt,nyoukado,nyougom,nyouj,nsaika
integer,allocatable :: nsety(:)
allocate(nsety(0:ny1*nsou+1))
!---ポアソン比
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poizx*eyy/ezz
poiyz=poiyx
!---せん断弾性係数
gxy=5.0d3/15
gxz=gxy
gyz=gxy
!---角度(1/tanθ)
kaku=(b4-b3)/2/(h1*7.d0/8)
!---要素長さ
dx1=(b1-b4)/2/nx1
dx2=b2/nx2
nx3=ny1*nsou*(nsou-1)/nsou!7h/8の高さに鋼材
dx3=(b4-b3)/2/nx3
dx4=b3/nx4
dy1=h1/(ny1*nsou)
dy2=h2/ny2
dz1=ell1/nz1
dz2=ell2/nz2
dz3=ell3/nz3
nz=nz1+nz2+nz3
!---節点座標の計算
nset=0
nsety=0
do k=0,nz1!内側
   do l=0,nsou-1
      if (l==0) then
         j1=0
      else
         j1=1
      endif
      do j=j1,ny1
         do i=0,nx2!H形鋼からはみ出たところ1
            nset=1+nset
            nsety(l*ny1+j+1)=nsety(l*ny1+j+1)+1
            x(nset)=i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            nsety(l*ny1+j+1)=nsety(l*ny1+j+1)+1
            x(nset)=b2+i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
         if (l==0) then
            i2=nx3!1番上の層は斜めにならない
         else
            i2=nint(nx3*(1.d0-real(l*ny1+j-ny1)/(ny1*(nsou-1))))!上から2~8層が斜めになる
         endif
         do i=1,i2!斜めのところ1
            nset=1+nset
            nsety(l*ny1+j+1)=nsety(l*ny1+j+1)+1
            x(nset)=b2+(b1-b4)/2+i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
         do i=0,nx2!H形鋼からはみ出たところ2
            nset=1+nset
            x(nset)=b1+b2*2.d0-i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            x(nset)=b2+b1-i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
         do i=1,i2!斜めのところ2
            nset=1+nset
            x(nset)=b2+(b1+b4)/2-i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=k*dz1
         enddo
      enddo
   enddo
enddo
nsety=nsety/(nz1+1)
do k=1,nz2!H形鋼の内側
   do l=0,nsou-1
      if (l==0) then
         j1=0
      else
         j1=1
      endif
      do j=j1,ny1
         do i=0,nx2!H形鋼からはみ出たところ1
            nset=1+nset
            x(nset)=i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            x(nset)=b2+i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
         if (l==0) then
            i2=nx3!1番上の層は斜めにならない
         else
            i2=nint(nx3*(1.d0-real(l*ny1+j-ny1)/(ny1*(nsou-1))))!上から2~8層が斜めになる
         endif
         do i=1,i2!斜めのところ1
            nset=1+nset
            x(nset)=b2+(b1-b4)/2+i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
         do i=0,nx2!H形鋼からはみ出たところ2
            nset=1+nset
            x(nset)=b1+b2*2.d0-i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            x(nset)=b2+b1-i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
         do i=1,i2!斜めのところ2
            nset=1+nset
            x(nset)=b2+(b1+b4)/2-i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+k*dz2
         enddo
      enddo
   enddo
enddo
do k=1,nz3!H形鋼の外側
   do l=0,nsou-1
      if (l==0) then
         j1=0
      else
         j1=1
      endif
      do j=j1,ny1
         do i=0,nx2!H形鋼からはみ出たところ1
            nset=1+nset
            x(nset)=i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            x(nset)=b2+i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
         if (l==0) then
            i2=nx3!1番上の層は斜めにならない
         else
            i2=nint(nx3*(1.d0-real(l*ny1+j-ny1)/(ny1*(nsou-1))))!上から2~8層が斜めになる
         endif
         do i=1,i2!斜めのところ1
            nset=1+nset
            x(nset)=b2+(b1-b4)/2+i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
         do i=0,nx2!H形鋼からはみ出たところ2
            nset=1+nset
            x(nset)=b1+b2*2.d0-i*dx2
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
         do i=1,nx1!H形鋼があるところ2
            nset=1+nset
            x(nset)=b2+b1-i*dx1
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
         do i=1,i2!斜めのところ2
            nset=1+nset
            x(nset)=b2+(b1+b4)/2-i*dx3
            y(nset)=l*h1/nsou+j*dy1
            z(nset)=ell1+ell2+k*dz3
         enddo
      enddo
   enddo
enddo
nsetclt=nset
!---ゴムの節点座標
do k=0,nz2!H形鋼の内側
   do j=1,ny2
      do i=0,nx1!H形鋼があるところ1
         nset=1+nset
         x(nset)=b2+i*dx1
         y(nset)=h1+j*dy2
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx3!斜めのところ1
         nset=1+nset
         x(nset)=b2+(b1-b4)/2+i*dx3
         y(nset)=h1+j*dy2
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx4!間
         nset=1+nset
         x(nset)=b2+(b1-b3)/2+i*dx4
         y(nset)=h1+j*dy2
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx3!斜めのところ2
         nset=1+nset
         x(nset)=b2+(b1+b3)/2+i*dx3
         y(nset)=h1+j*dy2
         z(nset)=ell1+k*dz2
      enddo
      do i=1,nx1!H形鋼があるところ2
         nset=1+nset
         x(nset)=b2+(b1+b4)/2+i*dx1
         y(nset)=h1+j*dy2
         z(nset)=ell1+k*dz2
      enddo
   enddo
enddo
do k=1,nz3!H形鋼の外側
   do j=1,ny2
      do i=0,nx1!H形鋼があるところ1
         nset=1+nset
         x(nset)=b2+i*dx1
         y(nset)=h1+j*dy2
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx3!斜めのところ1
         nset=1+nset
         x(nset)=b2+(b1-b4)/2+i*dx3
         y(nset)=h1+j*dy2
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx4!間
         nset=1+nset
         x(nset)=b2+(b1-b3)/2+i*dx4
         y(nset)=h1+j*dy2
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx3!斜めのところ2
         nset=1+nset
         x(nset)=b2+(b1+b3)/2+i*dx3
         y(nset)=h1+j*dy2
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx1!H形鋼があるところ2
         nset=1+nset
         x(nset)=b2+(b1+b4)/2+i*dx1
         y(nset)=h1+j*dy2
         z(nset)=ell1+ell2+k*dz3
      enddo
   enddo
enddo
nsetgom=nset-nsetclt
!---接続金物の節点座標
do k=0,nz1!内側
   j=ny1!上側
   do i=1,nx4-1
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=k*dz1
   enddo
   j=nsou*ny1!下側
   do i=1,nx3!斜めのところ１
      nset=1+nset
      x(nset)=b2+(b1-b4)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=k*dz1
   enddo
   do i=1,nx4!間
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=k*dz1
   enddo
   do i=1,nx3-1!斜めのところ2
      nset=1+nset
      x(nset)=b2+(b1+b3)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=k*dz1
   enddo
enddo
do k=1,nz2!H形鋼の内側
   j=ny1!上側
   do i=1,nx4-1
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=ell1+k*dz2
   enddo
   j=nsou*ny1!下側
   do i=1,nx3!斜めのところ１
      nset=1+nset
      x(nset)=b2+(b1-b4)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=ell1+k*dz2
   enddo
   do i=1,nx4!間
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=ell1+k*dz2
   enddo
   do i=1,nx3-1!斜めのところ2
      nset=1+nset
      x(nset)=b2+(b1+b3)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=ell1+k*dz2
   enddo
enddo
do k=1,nz3!H形鋼の外側
   j=ny1!上側
   do i=1,nx4-1
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=ell1+ell2+k*dz3
   enddo
   j=nsou*ny1!下側
   do i=1,nx3!斜めのところ１
      nset=1+nset
      x(nset)=b2+(b1-b4)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=ell1+ell2+k*dz3
   enddo
   do i=1,nx4!間
      nset=1+nset
      x(nset)=b2+(b1-b3)/2+i*dx4
      y(nset)=j*dy1
      z(nset)=ell1+ell2+k*dz3
   enddo
   do i=1,nx3-1!斜めのところ2
      nset=1+nset
      x(nset)=b2+(b1+b3)/2+i*dx3
      y(nset)=j*dy1
      z(nset)=ell1+ell2+k*dz3
   enddo
enddo
!---要素を構成する節点の計算
nyou=0
do k=0,nz-1!CLT
   do l=0,nsou-1
      do j=1,ny1
         do i=1,nsety(l*ny1+j+1)-1!左
            nyou=1+nyou
            mesh(nyou,1)=i+sum(nsety(0:l*ny1+j-1))*2+k*sum(nsety)*2
            mesh(nyou,2)=mesh(nyou,1)+nsety(l*ny1+j)*2
            mesh(nyou,3)=mesh(nyou,2)+sum(nsety)*2
            mesh(nyou,4)=mesh(nyou,1)+sum(nsety)*2
            mesh(nyou,5)=mesh(nyou,1)+1
            mesh(nyou,6)=mesh(nyou,2)+1
            mesh(nyou,7)=mesh(nyou,3)+1
            mesh(nyou,8)=mesh(nyou,4)+1
         enddo
         do i=1,nsety(l*ny1+j+1)-1!右
            nyou=1+nyou
            mesh(nyou,1)=nsety(l*ny1+j)+1+i+sum(nsety(0:l*ny1+j-1))*2+k*sum(nsety)*2
            mesh(nyou,2)=mesh(nyou,1)+nsety(l*ny1+j)+nsety(l*ny1+j+1)
            mesh(nyou,3)=mesh(nyou,2)+sum(nsety)*2
            mesh(nyou,4)=mesh(nyou,1)+sum(nsety)*2
            mesh(nyou,5)=mesh(nyou,1)-1
            mesh(nyou,6)=mesh(nyou,2)-1
            mesh(nyou,7)=mesh(nyou,3)-1
            mesh(nyou,8)=mesh(nyou,4)-1
         enddo
      enddo
   enddo
enddo
nyouclt=nyou
nyoukado=0
do k=0,nz-1!CLT角
   do l=1,nsou-1
      do j=1,ny1
         i=nsety(l*ny1+j+1)
         nyoukado=1+nyoukado!左
         meshk(nyoukado,1)=i+sum(nsety(0:l*ny1+j-1))*2+k*sum(nsety)*2
         meshk(nyoukado,2)=meshk(nyoukado,1)+1
         meshk(nyoukado,3)=meshk(nyoukado,1)+nsety(l*ny1+j)*2
         meshk(nyoukado,4)=meshk(nyoukado,1)+sum(nsety)*2
         meshk(nyoukado,5)=meshk(nyoukado,2)+sum(nsety)*2
         meshk(nyoukado,6)=meshk(nyoukado,3)+sum(nsety)*2
         nyoukado=1+nyoukado!右
         meshk(nyoukado,1)=nsety(l*ny1+j)+i+sum(nsety(0:l*ny1+j-1))*2+k*sum(nsety)*2
         meshk(nyoukado,2)=meshk(nyoukado,1)+nsety(l*ny1+j)+nsety(l*ny1+j+1)
         meshk(nyoukado,3)=meshk(nyoukado,1)+1
         meshk(nyoukado,4)=meshk(nyoukado,1)+sum(nsety)*2
         meshk(nyoukado,5)=meshk(nyoukado,2)+sum(nsety)*2
         meshk(nyoukado,6)=meshk(nyoukado,3)+sum(nsety)*2
      enddo
   enddo
enddo
nyougom=0
do k=0,nz2+nz3-1!ゴム
   do i=1,nx1!CLTとの境目の要素1
      nyougom=1+nyougom
      meshg(nyougom,1)=((sum(nsety)*2)*(nz1+k))+(sum(nsety(1:ny1*nsou))*2)+(nx2+i)
      meshg(nyougom,2)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+i
      meshg(nyougom,3)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*(1+k)+i
      meshg(nyougom,4)=meshg(nyougom,1)+sum(nsety)*2
      meshg(nyougom,5)=meshg(nyougom,1)+1
      meshg(nyougom,6)=meshg(nyougom,2)+1
      meshg(nyougom,7)=meshg(nyougom,3)+1
      meshg(nyougom,8)=meshg(nyougom,4)+1
   enddo
   do i=1,nx1!CLTとの境目の要素2
      nyougom=1+nyougom
      meshg(nyougom,1)=((sum(nsety)*2)*(nz1+k))+(sum(nsety(1:ny1*nsou))*2)+(nx2*2+nx1+2+i)
      meshg(nyougom,2)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1*2+nx3*2+nx4+1)-i
      meshg(nyougom,3)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*(k+1)+(nx1*2+nx3*2+nx4+1)-i
      meshg(nyougom,4)=meshg(nyougom,1)+sum(nsety)*2
      meshg(nyougom,5)=meshg(nyougom,1)-1
      meshg(nyougom,6)=meshg(nyougom,2)+1
      meshg(nyougom,7)=meshg(nyougom,3)+1
      meshg(nyougom,8)=meshg(nyougom,4)-1
   enddo
   do j=0,ny2-2!ゴムだけの要素1
      do i=1,nx1*2+nx3*2+nx4
         nyougom=1+nyougom
         !meshg(nyougom,1)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1+1)*j+i
         meshg(nyougom,1)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1*2+nx3*2+nx4+1)*j+i
         meshg(nyougom,2)=meshg(nyougom,1)+(nx1*2+nx3*2+nx4+1)
         meshg(nyougom,3)=meshg(nyougom,2)+(nx1*2+nx3*2+nx4+1)*ny2
         meshg(nyougom,4)=meshg(nyougom,1)+(nx1*2+nx3*2+nx4+1)*ny2
         meshg(nyougom,5)=meshg(nyougom,1)+1
         meshg(nyougom,6)=meshg(nyougom,2)+1
         meshg(nyougom,7)=meshg(nyougom,3)+1
         meshg(nyougom,8)=meshg(nyougom,4)+1
      enddo
   enddo
   do i=2,nx3*2+nx4-1!ゴムだけの要素2
      nyougom=1+nyougom
      meshg(nyougom,1)=nsetclt+nsetgom+(nx3*2+nx4*2-2)*(k+nz1)+(nx4-2+i)
      meshg(nyougom,2)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1+i)
      meshg(nyougom,3)=meshg(nyougom,2)+(nx1*2+nx3*2+nx4+1)*ny2
      meshg(nyougom,4)=meshg(nyougom,1)+(nx3*2+nx4*2-2)
      meshg(nyougom,5)=meshg(nyougom,1)+1
      meshg(nyougom,6)=meshg(nyougom,2)+1
      meshg(nyougom,7)=meshg(nyougom,3)+1
      meshg(nyougom,8)=meshg(nyougom,4)+1
   enddo
   nyougom=1+nyougom!CLTと接続金物との境界の要素1
   meshg(nyougom,1)=((sum(nsety)*2)*(nz1+k))+(sum(nsety(1:ny1*nsou))*2)+(nx2+nx1+1)
   meshg(nyougom,2)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1+1)
   meshg(nyougom,3)=meshg(nyougom,2)+((nx1*2+nx3*2+nx4+1)*ny2)
   meshg(nyougom,4)=meshg(nyougom,1)+sum(nsety)*2
   meshg(nyougom,5)=nsetclt+nsetgom+(nx3*2+nx4*2-2)*(k+nz1)+nx4
   meshg(nyougom,6)=meshg(nyougom,2)+1
   meshg(nyougom,7)=meshg(nyougom,3)+1
   meshg(nyougom,8)=meshg(nyougom,5)+(nx3*2+nx4*2-2)
   nyougom=1+nyougom!CLTと接続金物との境界の要素2
   meshg(nyougom,1)=nsetclt+nsetgom+(nx3*2+nx4*2-2)*(k+nz1)+(nx4*2+nx3*2-2)
   meshg(nyougom,2)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1+nx3*2+nx4)
   meshg(nyougom,3)=meshg(nyougom,2)+(nx1*2+nx3*2+nx4+1)*ny2
   meshg(nyougom,4)=meshg(nyougom,1)+(nx3*2+nx4*2-2)
   meshg(nyougom,5)=((sum(nsety)*2)*(nz1+k+1))
   meshg(nyougom,6)=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1+nx3*2+nx4+1)
   meshg(nyougom,7)=meshg(nyougom,6)+(nx1*2+nx3*2+nx4+1)*ny2
   meshg(nyougom,8)=meshg(nyougom,5)+sum(nsety)*2
enddo
nyouj=0!接続金物
do k=0,(nz1+nz2+nz3)/2-1
   nyouj=1+nyouj!CLTとの間1
   meshj(nyouj,1)=sum(nsety(1:ny1))*2+nsety((ny1-1)+2)+(sum(nsety)*2*2)*k
   meshj(nyouj,2)=nsetclt+nsetgom+((nx4*2+nx3*2-2)*2)*k+2
   meshj(nyouj,3)=meshj(nyouj,2)+(nx4*2+nx3*2-2)*2
   meshj(nyouj,4)=meshj(nyouj,1)+(sum(nsety)*2*2)
   meshj(nyouj,5)=meshj(nyouj,2)-1
   meshj(nyouj,6)=meshj(nyouj,2)+(nx4*2+nx3*2-2)
   meshj(nyouj,7)=meshj(nyouj,3)-1
   meshj(nyouj,8)=meshj(nyouj,1)+(sum(nsety)*2)
   do i=1,(nx4-4)/2!鋼材だけの要素1
      nyouj=1+nyouj
      meshj(nyouj,1)=nsetclt+nsetgom+((nx4*2+nx3*2-2)*2)*k+2*i
      meshj(nyouj,2)=meshj(nyouj,1)+2
      meshj(nyouj,3)=meshj(nyouj,2)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,4)=meshj(nyouj,1)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,5)=meshj(nyouj,1)+1
      meshj(nyouj,6)=meshj(nyouj,2)+(nx4*2+nx3*2-2)
      meshj(nyouj,7)=meshj(nyouj,5)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,8)=meshj(nyouj,1)+(nx4*2+nx3*2-2)
   enddo
   nyouj=1+nyouj!CLTとの間2
   meshj(nyouj,1)=nsetclt+nsetgom+((nx4*2+nx3*2-2)*2)*k+(nx4-2)
   meshj(nyouj,2)=sum(nsety(1:ny1))*2+nsety((ny1-1)+2)*2+(sum(nsety)*2*2)*k
   meshj(nyouj,3)=meshj(nyouj,2)+(sum(nsety)*2*2)
   meshj(nyouj,4)=meshj(nyouj,1)+(nx4*2+nx3*2-2)*2
   meshj(nyouj,5)=meshj(nyouj,1)+1
   meshj(nyouj,6)=meshj(nyouj,2)+(sum(nsety)*2)
   meshj(nyouj,7)=meshj(nyouj,4)+1
   meshj(nyouj,8)=meshj(nyouj,1)+(nx4*2+nx3*2-2)
   nyouj=1+nyouj!CLTとの間3
   meshj(nyouj,1)=sum(nsety(1:ny1+nx3))*2+nsety(ny1+nx3+1)+(sum(nsety)*2*2)*k
   meshj(nyouj,2)=nsetclt+nsetgom+(nx4-1)+((nx4*2+nx3*2-2)*2)*k+2
   meshj(nyouj,3)=meshj(nyouj,2)+(nx4*2+nx3*2-2)*2
   meshj(nyouj,4)=meshj(nyouj,1)+(sum(nsety)*2*2)
   meshj(nyouj,5)=meshj(nyouj,2)-1
   meshj(nyouj,6)=meshj(nyouj,2)+(nx4*2+nx3*2-2)
   meshj(nyouj,7)=meshj(nyouj,3)-1
   meshj(nyouj,8)=meshj(nyouj,1)+(sum(nsety)*2)
   do i=1,(nx3*2+nx4-4)/2!鋼材だけの要素2
      nyouj=1+nyouj
      meshj(nyouj,1)=nsetclt+nsetgom+(nx4-1)+((nx4*2+nx3*2-2)*2)*k+2*i
      meshj(nyouj,2)=meshj(nyouj,1)+2
      meshj(nyouj,3)=meshj(nyouj,2)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,4)=meshj(nyouj,1)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,5)=meshj(nyouj,1)+1
      meshj(nyouj,6)=meshj(nyouj,2)+(nx4*2+nx3*2-2)
      meshj(nyouj,7)=meshj(nyouj,5)+(nx4*2+nx3*2-2)*2
      meshj(nyouj,8)=meshj(nyouj,1)+(nx4*2+nx3*2-2)
   enddo
   nyouj=1+nyouj!CLTとの間4
   meshj(nyouj,1)=nsetclt+nsetgom+(nx4*2+nx3*2-2)+((nx4*2+nx3*2-2)*2)*k-1
   meshj(nyouj,2)=sum(nsety)*2+(sum(nsety)*2*2)*k
   meshj(nyouj,3)=meshj(nyouj,2)+(sum(nsety)*2*2)
   meshj(nyouj,4)=meshj(nyouj,1)+(nx4*2+nx3*2-2)*2
   meshj(nyouj,5)=meshj(nyouj,1)+1
   meshj(nyouj,6)=meshj(nyouj,2)+(sum(nsety)*2)
   meshj(nyouj,7)=meshj(nyouj,4)+1
   meshj(nyouj,8)=meshj(nyouj,1)+(nx4*2+nx3*2-2)
   do i=1,nx3/2!CLTと重なっている要素1
      nyouj=1+nyouj
      meshj(nyouj,1)=sum(nsety(1:ny1+2*(i-1)))*2+nsety((ny1-1)+2*i)+(sum(nsety)*2*2)*k
      meshj(nyouj,2)=meshj(nyouj,1)+(nsety((ny1-1)+2*i)+nsety(ny1+2*i)*2+nsety((ny1+1)+2*i))
      meshj(nyouj,3)=meshj(nyouj,2)+(sum(nsety)*2*2)
      meshj(nyouj,4)=meshj(nyouj,1)+(sum(nsety)*2*2)
      meshj(nyouj,5)=meshj(nyouj,1)+(nsety((ny1-1)+2*i)+nsety(ny1+2*i))
      meshj(nyouj,6)=meshj(nyouj,2)+(sum(nsety)*2)
      meshj(nyouj,7)=meshj(nyouj,5)+(sum(nsety)*2*2)
      meshj(nyouj,8)=meshj(nyouj,1)+(sum(nsety)*2)
   enddo
   do i=1,nx3/2!CLTと重なっている要素2
      nyouj=1+nyouj
      meshj(nyouj,1)=sum(nsety(1:ny1+2*(i-1)))*2+nsety((ny1-1)+2*i)*2+(sum(nsety)*2*2)*k
      meshj(nyouj,2)=meshj(nyouj,1)+(sum(nsety)*2*2)
      meshj(nyouj,3)=meshj(nyouj,2)+(nsety(ny1+2*i)*2+nsety((ny1+1)+2*i)*2)
      meshj(nyouj,4)=meshj(nyouj,1)+(nsety(ny1+2*i)*2+nsety((ny1+1)+2*i)*2)
      meshj(nyouj,5)=meshj(nyouj,1)+(sum(nsety)*2)
      meshj(nyouj,6)=meshj(nyouj,2)+(nsety(ny1+2*i)*2)
      meshj(nyouj,7)=meshj(nyouj,4)+(sum(nsety)*2)
      meshj(nyouj,8)=meshj(nyouj,1)+(nsety(ny1+2*i)*2)
   enddo

enddo
!---節点座標の出力
print'("*NODE, NSET=Nall")'
do i=1,nset
   print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
       &i,',',x(i),',',y(i),',',z(i)
end do
!---要素の出力
do k=0,nsou-1!CLT
   print'(A,i1)',"*ELEMENT, TYPE=C3D8, ELSET=E",k+1
   do j=0,nz-1
      do i=1,(sum(nsety(k*ny1+2:(k+1)*ny1+1))-ny1)*2
         l=i+j*(sum(nsety)-nsety(1)-ny1*nsou)*2+&
             &((sum(nsety(1:k*ny1+1))-nsety(1))-(ny1*k))*2
         print'(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
do k=0,nsou-2!CLT角
   print'(A,i1)',"*ELEMENT, TYPE=C3D6, ELSET=Ekado",k+2
   do j=0,nz-1
      do i=1,ny1*2
         l=i+j*(ny1*2*(nsou-1))+k*(ny1*2)
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
             &nyouclt+l,',',&
             &meshk(l,1),',',&
             &meshk(l,2),',',&
             &meshk(l,3),',',&
             &meshk(l,4),',',&
             &meshk(l,5),',',&
             &meshk(l,6)
      enddo
   enddo
enddo
print'(A)',"*ELEMENT, TYPE=C3D8, ELSET=Egom"!ゴム
do l=1,nyougom
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
       &l+nyouclt+nyoukado,',',&
       &meshg(l,1),',',&
       &meshg(l,2),',',&
       &meshg(l,3),',',&
       &meshg(l,4),',',&
       &meshg(l,5),',',&
       &meshg(l,6),',',&
       &meshg(l,7),',',&
       &meshg(l,8)
enddo
print'(A)',"*ELEMENT, TYPE=S8, ELSET=Ekou1"!接続金物1
do j=0,nz/2-1
   do i=1,nx4/2
      l=(nx4+nx3*2)*j+i
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
       &l+nyouclt+nyoukado+nyougom,',',&
       &meshj(l,1),',',&
       &meshj(l,2),',',&
       &meshj(l,3),',',&
       &meshj(l,4),',',&
       &meshj(l,5),',',&
       &meshj(l,6),',',&
       &meshj(l,7),',',&
       &meshj(l,8)
   enddo
enddo
print'(A)',"*ELEMENT, TYPE=S8, ELSET=Ekou2"!接続金物2
do j=0,nz/2-1
   do i=1,nx3+nx4/2
      l=nx4/2+(nx4+nx3*2)*j+i
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
       &l+nyouclt+nyoukado+nyougom,',',&
       &meshj(l,1),',',&
       &meshj(l,2),',',&
       &meshj(l,3),',',&
       &meshj(l,4),',',&
       &meshj(l,5),',',&
       &meshj(l,6),',',&
       &meshj(l,7),',',&
       &meshj(l,8)
   enddo
enddo
print'(A)',"*ELEMENT, TYPE=S8, ELSET=Ekou3"!接続金物3
do j=0,nz/2-1
   do i=1,nx3
      l=nx4+nx3+(nx4+nx3*2)*j+i
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
       &l+nyouclt+nyoukado+nyougom,',',&
       &meshj(l,1),',',&
       &meshj(l,2),',',&
       &meshj(l,3),',',&
       &meshj(l,4),',',&
       &meshj(l,5),',',&
       &meshj(l,6),',',&
       &meshj(l,7),',',&
       &meshj(l,8)
   enddo
enddo
!#--境界条件の節点群
print'("*NSET,NSET=kotei1")'!H形鋼の全面
do k=0,nz2+nz3
   do i=1,nx1*2+nx3*2+nx4+1
      l=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1*2+nx3*2+nx4+1)+i
      print'(i6,a)',l
   end do
enddo
print'("*NSET,NSET=kotei2")'!H形鋼の半面
do k=0,nz2
   do i=1,nx1*2+nx3*2+nx4+1
      l=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1*2+nx3*2+nx4+1)+i
      print'(i6,a)',l
   end do
enddo
print'("*NSET,NSET=kotei3")'!ボルトの線上
k=nz2
do i=1,nx1*2+nx3*2+nx4+1
   l=nsetclt+((nx1*2+nx3*2+nx4+1)*ny2)*k+(nx1*2+nx3*2+nx4+1)+i
   print'(i6,a)',l
end do
print'("*NSET,NSET=koteiz")'!z軸対称面
do i=1,nsetclt
   if (z(i)>-dz1/2.and.z(i)<dz1/2) then
      print'(i6,a)',i
   endif
end do
!---載荷条件の節点群
print'("*NSET,NSET=saika")'
nsaika=0
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>itiz-saikaz/2-dz1/3.and.z(i)<itiz+saikaz/2+dz1/3.and.&
       &x(i)>itix-saikax/2-dx3/3.and.x(i)<itix+saikax/2+dx3/3) then
      nsaika=nsaika+1
      print'(i6,a)',i
   endif
end do
!#--変位計の節点群
print'("*NSET,NSET=hen1")'!変位計1
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>250.d-3-dx2/2.and.x(i)<250.d-3+dx2/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen2")'!変位計2
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>1000.d-3-dx1/2.and.x(i)<1000.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen3")'!変位計3
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>1750.d-3-dx1/2.and.x(i)<1750.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen7")'!変位計7
do i=nsetclt+nsetgom+1,nset
   if (y(i)>h1/nsou-dy1/3.and.y(i)<h1/nsou+dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>2010.d-3-dx4/2.and.x(i)<2010.d-3+dx4/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen4")'!変位計4
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>2250.d-3-dx1/2.and.x(i)<2250.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen5")'!変位計5
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>3000.d-3-dx1/2.and.x(i)<3000.d-3+dx1/2) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=hen6")'!変位計6
do i=1,nsetclt
   if (y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>-dz1/2.and.z(i)<dz1/2.and.&
       &x(i)>3750.d-3-dx2/2.and.x(i)<3750.d-3+dx2/2) then
      print'(i6,a)',i
   endif
end do
!#--境界条件
print'("*BOUNDARY")'
print'("kotei3,1,3")'
print'("koteiz,3")'!z軸対称面
!---材料設定
print'("*MATERIAL,NAME=EL3")'!鋼材
print'("*ELASTIC")'!等方性弾性体
print'(f12.4,a,f12.4)',ekou,',',poikou
print'("*MATERIAL,NAME=EL2")'!ゴム
print'("*ELASTIC")'!等方性弾性体
print'(f12.4,a,f12.4)',egom,',',poigom
print'("*MATERIAL,NAME=EL1")'!CLT
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'("*DENSITY")'!密度
!print'("0.38E-9")'
!#--局所座標系の設定
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'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E2,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E3,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E4,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E5,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E6,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E7,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E8,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado2,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado3,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado4,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado5,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado6,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado7,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ekado8,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Egom,MATERIAL=EL2")'!ゴム
print'("*SHELL SECTION,ELSET=Ekou1,MATERIAL=EL3,OFFSET=-0.009")'!接続金物(上),オフセット量を調整すべし！
print'(1pe13.5)',t
print'("*SHELL SECTION,ELSET=Ekou2,MATERIAL=EL3,OFFSET=0.003")'!接続金物(下),オフセット量を調整すべし！
print'(1pe13.5)',t
print'("*SHELL SECTION,ELSET=Ekou3,MATERIAL=EL3,OFFSET=0.008")'!接続金物(斜),オフセット量を調整すべし！
print'(1pe13.5)',t
!---載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nsaika
!print'("*FREQUENCY")'!固有値解析
!print'("10")'!振動モード
!#--datファイルの設定
print'("*NODE PRINT,NSET=hen1")'!変位計1
print'("U")'
print'("*NODE PRINT,NSET=hen2")'!変位計2
print'("U")'
print'("*NODE PRINT,NSET=hen3")'!変位計3
print'("U")'
print'("*NODE PRINT,NSET=hen7")'!変位計7
print'("U")'
print'("*NODE PRINT,NSET=hen4")'!変位計4
print'("U")'
print'("*NODE PRINT,NSET=hen5")'!変位計5
print'("U")'
print'("*NODE PRINT,NSET=hen6")'!変位計6
print'("U")'
print'("*NODE PRINT,NSET=saika")'!平均化したほうがいい？
print'("U")'
!print'("*EL PRINT,ELSET=")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
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
