!輪荷重を受けたCLT床版の載荷試験(CCX)
! 1/4解析
! ボルトを梁要素として混在させる。
! 横桁のH鋼は変形しないものと仮定してモデリングせず。
implicit none
!#--入力
real    :: b1=2010.d-3    ! m   :支間方向の幅(x)
real    :: h=240.d-3      ! m   :全層の高さ(y)
real    :: ell1=400.d-3   ! m   :橋軸方向の長さ1(z)
real    :: ell2=1350.d-3  ! m   :橋軸方向の長さ2(z)
real    :: ell3=500.d-3   ! m   :橋軸方向の長さ3(z)
real    :: kotei=80.d-3   ! m   :端からボルトまでの距離(x)
real    :: d=16.d-3       ! m   :ボルトの直径
integer :: nb=4           ! m   :ボルトの本数(1/4)
real    :: b2=300.d-3     ! m   :床版を固定しているH型鋼の幅
real    :: ellh=3380.d-3  ! m   :床版を固定しているH型鋼の長さ
integer :: nsou=8         !     :層数
integer :: nx1=2!4          !     :支間方向分割数1
integer :: nx2=80!120        !     :支間方向分割数2
integer :: ny=4!4           !偶数 :1層ごとの高さ方向分割数
integer :: nz1=12!18         !     :橋軸方向分割数1
integer :: nz2=20!32         !     :橋軸方向分割数2(1ボルト間)
integer :: nz3=12!16         !     :橋軸方向分割数3
real    :: p=9.58d-3      ! MN  :荷重
real    :: saikax=500.d-3 ! m   :載荷板の支間方向幅
real    :: saikaz=200.d-3 ! m   :載荷板の橋軸幅
real    :: itix=1005.d-3  ! m   :支間方向(x)の載荷板中心位置
real    :: itiz=2000.d-3  ! m   :橋軸方向(z)の載荷板中心位置
real    :: exx=5.0d3/25   ! MPa :x軸方向のヤング率
real    :: eyy=5.0d3/25   ! MPa :y軸方向のヤング率
real    :: ezz=5.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 :: x(1000000),y(1000000),z(1000000)     !節点座標
integer :: mesh1(1000000,8),mesh2(10000,3) !要素と要素を構成する節点
real :: e66(6,6),b(3,3),dx1,dx2,dy,dz1,dz2,dz3,ell
integer i,j,k,l,nset,nyou1,nyou2,nsaika,nx,nz
!---ポアソン比
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poixy
poiyz=poixy
!---せん断弾性係数
gxy=ezz/15
gxz=gxy
gyz=gxy
!---要素長さ
dx1=kotei/nx1
dx2=(b1-kotei*2.d0)/2/nx2
dy=h/(ny*nsou)
dz1=ell1/nz1
dz2=ell2/(nb-1)/nz2
dz3=ell3/2/nz3
!---要素分割数や長さ合計
nz=nz1+nz2*(nb-1)+nz3
nx=nx1+nx2
ell=ell1+ell2+ell3/2
!---節点座標の計算
nset=0
do k=0,nz1
   do j=0,ny*nsou
      do i=0,nx1
         nset=1+nset
         x(nset)=i*dx1
         y(nset)=j*dy
         z(nset)=k*dz1
      enddo
      do i=1,nx2
         nset=1+nset
         x(nset)=kotei+i*dx2
         y(nset)=j*dy
         z(nset)=k*dz1
      enddo
   enddo
enddo
do l=0,nb-2
   do k=1,nz2
      do j=0,ny*nsou
         do i=0,nx1
            nset=1+nset
            x(nset)=i*dx1
            y(nset)=j*dy
            z(nset)=ell1+ell2/(nb-1)*l+k*dz2
         enddo
         do i=1,nx2
            nset=1+nset
            x(nset)=kotei+i*dx2
            y(nset)=j*dy
            z(nset)=ell1+ell2/(nb-1)*l+k*dz2
         enddo
      enddo
   enddo
enddo
do k=1,nz3
   do j=0,ny*nsou
      do i=0,nx1
         nset=1+nset
         x(nset)=i*dx1
         y(nset)=j*dy
         z(nset)=ell1+ell2+k*dz3
      enddo
      do i=1,nx2
         nset=1+nset
         x(nset)=kotei+i*dx2
         y(nset)=j*dy
         z(nset)=ell1+ell2+k*dz3
      enddo
   enddo
enddo
!---要素を構成する節点の計算(木材)
nyou1=0
do k=0,nz-1
   do j=0,ny*nsou-1
      do i=1,nx
         nyou1=1+nyou1
         mesh1(nyou1,1)=i+j*(nx+1)+k*((nx+1)*(ny*nsou+1))
         mesh1(nyou1,2)=mesh1(nyou1,1)+(nx+1)
         mesh1(nyou1,3)=mesh1(nyou1,2)+(nx+1)*(ny*nsou+1)
         mesh1(nyou1,4)=mesh1(nyou1,1)+(nx+1)*(ny*nsou+1)
         mesh1(nyou1,5)=mesh1(nyou1,1)+1
         mesh1(nyou1,6)=mesh1(nyou1,5)+nx+1
         mesh1(nyou1,7)=mesh1(nyou1,6)+(nx+1)*(ny*nsou+1)
         mesh1(nyou1,8)=mesh1(nyou1,5)+(nx+1)*(ny*nsou+1)
      enddo
   enddo
enddo
!---要素を構成する節点の計算(ボルト)
nyou2=0
do j=0,(ny*nsou)/2-1
   nyou2=1+nyou2
   mesh2(nyou2,1)=(nx+1)*(ny*nsou+1)*nz1+(nx1+1)+(nx+1)*2*j
   mesh2(nyou2,2)=mesh2(nyou2,1)+(nx+1)
   mesh2(nyou2,3)=mesh2(nyou2,2)+(nx+1)
enddo
do k=1,(nb-1)
   do j=0,(ny*nsou)/2-1
      nyou2=1+nyou2
      mesh2(nyou2,1)=(nx+1)*(ny*nsou+1)*(nz1+nz2*k)+(nx1+1)+(nx+1)*2*j
      mesh2(nyou2,2)=mesh2(nyou2,1)+(nx+1)
      mesh2(nyou2,3)=mesh2(nyou2,2)+(nx+1)
   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
   print'(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)
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
             &l,',',&
             &mesh1(l,1),',',&
             &mesh1(l,2),',',&
             &mesh1(l,3),',',&
             &mesh1(l,4),',',&
             &mesh1(l,5),',',&
             &mesh1(l,6),',',&
             &mesh1(l,7),',',&
             &mesh1(l,8)
      enddo
   enddo
enddo
!---要素の出力(ボルト)
print'(A,i1)',"*ELEMENT, TYPE=B32, ELSET=Eb"
do l=1,nyou2
   print'(i6,a,i6,a,i6,a,i6)',nyou1+l,',',mesh2(l,1),',',mesh2(l,2),',',mesh2(l,3)
enddo
!#--境界条件の節点群
print'("*NSET,NSET=kotei1")'!H鋼
do i=1,nset
   if (x(i)>kotei+d/2.and.x(i)<b2/2+dx2/3.and.&!H鋼の面半分で固定する場合
   !if (x(i)>-dx1/3.and.x(i)<b2/2+dx2/3.and.&!H鋼の全面で固定する場合
   !if (x(i)>kotei-dx1/3.and.x(i)<kotei+dx2/3.and.&!ボルトらへんで線で固定する場合
       &y(i)>h-dy/3.and.y(i)<h+dy/3.and.&
       &z(i)>ell-ellh/2-dz1/3.and.z(i)<ell+dz1/3) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=kotei2")'!ボルト
do i=ny*nsou/2,(ny*nsou/2)*nb,ny*nsou/2
   print'(i6,a)',mesh2(i,3)
enddo
print'("*NSET,NSET=kotei3")'!x対称面
do i=1,nset
   if (x(i)>b1/2-dx2/3.and.x(i)<b1/2+dx2/3) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=kotei4")'!z対称面
do i=1,nset
   if (z(i)>ell-dz3/3.and.z(i)<ell+dz3/3) then
      print'(i6,a)',i
   endif
end do
!---載荷条件の節点群
print'("*NSET,NSET=saika")'
nsaika=0
do i=1,nset
   if (y(i)>-dy/3.and.y(i)<dy/3.and.&
       &z(i)>=itiz-saikaz/2-dz3/3.and.z(i)<=itiz+dz3/3.and.&
       &x(i)>=itix-saikax/2-dx2/3.and.x(i)<=itix+dx2/3) then
      nsaika=nsaika+1
      print'(i6,a)',i
   endif
end do
!!#--ポイントゲージの節点群
!print'("*NSET,NSET=henni")'
!do i=1,nset
!   !---高さが3/7*hの場合
!   if (y(i)>h*3.d0/7-dy/2.and.y(i)<h*3.d0/7+dy/2.and.&
!       &z(i)>=ell/2-saikaz/2-dz/2.and.z(i)<=ell/2+saikaz/2+dz/2.and.&
!       &x(i)>=-dx2/3.and.x(i)<=dx2/3) then
!      print'(i6,a)',i
!   else if (y(i)>h*3.d0/7-dy/2.and.y(i)<h*3.d0/7+dy/2.and.&
!       &z(i)>=ell/2-saikaz/2-dz/2.and.z(i)<=ell/2+saikaz/2+dz/2.and.&
!       &x(i)>=b1-dx2/3.and.x(i)<=b1+dx2/3) then
!      print'(i6,a)',i
!   endif
!end do
!#--境界条件
print'("*BOUNDARY")'
print'("kotei1,2")'!H鋼(鉛直方向だけ固定)
!print'("kotei1,1,3")'!H鋼(全方向固定)
print'("kotei2,3")'!ボルト(橋軸方向固定:主桁は縮まないと思うので)
!print'("kotei2,1,3")'!ボルト(全方向固定:この場合*EQUATIONをコメントアウトする)
print'("kotei3,1")'!x対称面(滑る)
print'("kotei4,3")'!z対称面(滑る)
print'("*EQUATION")'!ボルトの支間方向変位を一様に等しく
do i=ny*nsou/2,(ny*nsou/2)*(nb-1),ny*nsou/2
   print'("2")'!H鋼1
   print'(i6,A,i6,A)',mesh2(i,3),',1,1,',mesh2(i+ny*nsou/2,3),',1,-1'
enddo
!---材料設定
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'("*DENSITY")'!密度
!print'("0.38E-9")'
print'("*MATERIAL,NAME=ELkou")'!鋼材
print'("*ELASTIC")'
print'(f12.4,a,f12.4)',ekou,',',poikou
!print'("*DENSITY")'!密度
!print'("7.85E-9")'
print'("*BEAM SECTION,ELSET=Eb,MATERIAL=ELkou,SECTION=CIRC")'
print'(f12.4,a,f12.4)',d,',',d
print'("0.,0.,1.")'
!#--局所座標系の設定
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=ELmoku,ORIENTATION=OR2")'
print'("*SOLID SECTION,ELSET=E2,MATERIAL=ELmoku,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E3,MATERIAL=ELmoku,ORIENTATION=OR2")'
print'("*SOLID SECTION,ELSET=E4,MATERIAL=ELmoku,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E5,MATERIAL=ELmoku,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E6,MATERIAL=ELmoku,ORIENTATION=OR2")'
print'("*SOLID SECTION,ELSET=E7,MATERIAL=ELmoku,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E8,MATERIAL=ELmoku,ORIENTATION=OR2")'
!---載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/4/nsaika
!print'("*FREQUENCY")'!固有値解析
!print'("10")'!振動モード
!#--datファイルの設定
print'("*NODE PRINT,NSET=saika")'
print'("U")'
!print'("*NODE PRINT,NSET=henni")'
!print'("U")'
print'("*EL PRINT,ELSET=Eb")'
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
