!中村先生のCLT橋の3点曲げ(CCX)
! CLTと集成材は節点共有なので，完全にずれない．
! ラグボルトも含んでいない．
implicit none
!---入力
!#--共通の寸法等
real    :: ell1=5800.d-3   ! m   :スパン(z)
real    :: ell2=100.d-3    ! m   :端から支点までの距離(z)
integer :: nz1=200         !     :軸方向分割数1
integer :: nz2=2           !     :軸方向分割数2
real    :: poizx=0.4d0     !     :ポアソン比
!#--床版(CLT)の寸法等
real    :: b1=1000.d-3     ! m   :幅(x)
real    :: h1=150.d-3      ! m   :全層の高さ(y)
integer :: nsou=5          !     :層数
integer :: nx1=8           !     :両脇の幅方向分割数
integer :: ny1=2           !     :1層ごとの高さ方向分割数
real    :: exx1=5.0d3/25    ! MPa :x軸(幅)方向のヤング率
real    :: eyy1=5.0d3/25    ! MPa :y軸(高さ)方向のヤング率
real    :: ezz1=5.0d3       ! MPa :z軸(長さ)方向のヤング率
!#--梁(集成材)の寸法等
real    :: b2=120.d-3      ! m   :幅(x)
real    :: h2=500.d-3      ! m   :高さ(y)
integer :: nx2=4           !     :梁の幅方向分割数
integer :: ny2=20          !     :高さ方向分割数
real    :: exx2=7.0d3/25    ! MPa :x軸(幅)方向のヤング率
real    :: eyy2=7.0d3/25    ! MPa :y軸(高さ)方向のヤング率
real    :: ezz2=7.0d3       ! MPa :z軸(長さ)方向のヤング率
!#--材料定数や境界条件，載荷条件など
real    :: p=150.d-3       ! MN  :荷重
real    :: saikax=1000.d-3 ! m   :載荷板の幅(x)
real    :: saikaz=50.d-3   ! m   :載荷板の幅(z)
real    :: itix=500.d-3    ! m   :幅方向(x)の載荷板中心位置
real    :: itiz=3000.d-3   ! m   :長さ方向(z)の載荷板中心位置
!---
real poixy,poixz,poiyx,poiyz,poizy,gxy1,gxz1,gyz1,gxy2,gxz2,gyz2
real :: x(1000000),y(1000000),z(1000000) !節点座標
integer :: mesh(1000000,8)               !要素と要素を構成する節点
real :: e66(6,6),b(3,3),dx1,dx2,dy1,dy2,dz1,dz2
integer i,j,k,l,nset,nset1,nset2,nyou,nsaika
!---ポアソン比
poizy=poizx
poixz=poizx*exx1/ezz1
poixy=poixz
poiyz=poizy*eyy1/ezz1
poiyx=poiyz
!---せん断弾性係数
gxy1=ezz1/15!CLT
gxz1=gxy1
gyz1=gxy1
gxy2=ezz2/15!集成材
gxz2=gxy2
gyz2=gxy2
!---要素長さ
dx1=(b1-b2)/2/nx1
dx2=b2/nx2
dy1=h1/(ny1*nsou)
dy2=h2/ny2
dz1=ell1/nz1
dz2=ell2/nz2
!---節点座標の計算
nset=0
do k=0,nz2!CLT(端1)
   do j=0,ny1*nsou
      do i=0,nx1!端1
         nset=1+nset
         x(nset)=i*dx1
         y(nset)=j*dy1
         z(nset)=k*dz2
      enddo
      do i=1,nx2!中
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=j*dy1
         z(nset)=k*dz2
      enddo
      do i=1,nx1!端2
         nset=1+nset
         x(nset)=nx1*dx1+nx2*dx2+i*dx1
         y(nset)=j*dy1
         z(nset)=k*dz2
      enddo
   enddo
enddo
do k=1,nz1!CLT(中)
   do j=0,ny1*nsou
      do i=0,nx1!端1
         nset=1+nset
         x(nset)=i*dx1
         y(nset)=j*dy1
         z(nset)=ell2+k*dz1
      enddo
      do i=1,nx2!中
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=j*dy1
         z(nset)=ell2+k*dz1
      enddo
      do i=1,nx1!端2
         nset=1+nset
         x(nset)=nx1*dx1+nx2*dx2+i*dx1
         y(nset)=j*dy1
         z(nset)=ell2+k*dz1
      enddo
   enddo
enddo
do k=1,nz2!CLT(端2)
   do j=0,ny1*nsou
      do i=0,nx1!端1
         nset=1+nset
         x(nset)=i*dx1
         y(nset)=j*dy1
         z(nset)=ell2+ell1+k*dz2
      enddo
      do i=1,nx2!中
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=j*dy1
         z(nset)=ell2+ell1+k*dz2
      enddo
      do i=1,nx1!端2
         nset=1+nset
         x(nset)=nx1*dx1+nx2*dx2+i*dx1
         y(nset)=j*dy1
         z(nset)=ell2+ell1+k*dz2
      enddo
   enddo
enddo
nset1=nset
do k=0,nz2!集成材(端1)
   do j=1,ny2
      do i=0,nx2
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=h1+j*dy2
         z(nset)=k*dz2
      enddo
   enddo
enddo
do k=1,nz1!集成材(中)
   do j=1,ny2
      do i=0,nx2
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=h1+j*dy2
         z(nset)=ell2+k*dz1
      enddo
   enddo
enddo
do k=1,nz2!集成材(端2)
   do j=1,ny2
      do i=0,nx2
         nset=1+nset
         x(nset)=nx1*dx1+i*dx2
         y(nset)=h1+j*dy2
         z(nset)=ell2+ell1+k*dz2
      enddo
   enddo
enddo
nset2=nset-nset1
!---要素を構成する節点の計算
nyou=0
do k=0,nz1+nz2*2-1!CLT
   do j=0,ny1*nsou-1
      do i=1,nx1*2+nx2
         nyou=1+nyou
         mesh(nyou,1)=i+j*((nx1*2+nx2)+1)+k*(((nx1*2+nx2)+1)*(ny1*nsou+1))
         mesh(nyou,2)=mesh(nyou,1)+((nx1*2+nx2)+1)
         mesh(nyou,3)=mesh(nyou,2)+((nx1*2+nx2)+1)*(ny1*nsou+1)
         mesh(nyou,4)=mesh(nyou,1)+((nx1*2+nx2)+1)*(ny1*nsou+1)
         mesh(nyou,5)=mesh(nyou,1)+1
         mesh(nyou,6)=mesh(nyou,5)+(nx1*2+nx2)+1
         mesh(nyou,7)=mesh(nyou,6)+((nx1*2+nx2)+1)*(ny1*nsou+1)
         mesh(nyou,8)=mesh(nyou,5)+((nx1*2+nx2)+1)*(ny1*nsou+1)
      enddo
   enddo
enddo
do k=0,nz1+nz2*2-1!集成材
   do i=1,nx2!CLTと節点を共有する要素
      nyou=1+nyou
      mesh(nyou,1)=((nx1*2+nx2+1)*(ny1*nsou)+nx1)+i+k*(((nx1*2+nx2)+1)*(ny1*nsou+1))
      mesh(nyou,2)=nset1+i+k*((nx2+1)*ny2)
      mesh(nyou,3)=mesh(nyou,2)+(nx2+1)*ny2
      mesh(nyou,4)=mesh(nyou,1)+((nx1*2+nx2)+1)*(ny1*nsou+1)
      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 j=0,ny2-2!集成材の節点のみの要素
      do i=1,nx2
         nyou=1+nyou
         mesh(nyou,1)=nset1+i+j*(nx2+1)+k*((nx2+1)*ny2)
         mesh(nyou,2)=mesh(nyou,1)+nx2+1
         mesh(nyou,3)=mesh(nyou,2)+(nx2+1)*ny2
         mesh(nyou,4)=mesh(nyou,1)+(nx2+1)*ny2
         mesh(nyou,5)=mesh(nyou,1)+1
         mesh(nyou,6)=mesh(nyou,5)+nx2+1
         mesh(nyou,7)=mesh(nyou,6)+(nx2+1)*ny2
         mesh(nyou,8)=mesh(nyou,5)+(nx2+1)*ny2
      enddo
   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,nz1+nz2*2-1
      do i=1,(nx1*2+nx2)*ny1
         l=i+j*((nx1*2+nx2)*ny1*nsou)+k*((nx1*2+nx2)*ny1)
         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
print'(A)',"*ELEMENT, TYPE=C3D8, ELSET=Ehari"!集成材
do i=l+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)
enddo
!---境界条件の節点群
print'("*NSET,NSET=kotei1")'
do i=nset1+1,nset
   if (y(i)>h1+h2-dy2/3.and.y(i)<h1+h2+dy2/3.and.&
       &z(i)>ell2-dz2/3.and.z(i)<ell2+dz2/3) then
      print'(i6,a)',i
   endif
end do
print'("*NSET,NSET=kotei2")'
do i=nset1+1,nset
   if (y(i)>h1+h2-dy2/3.and.y(i)<h1+h2+dy2/3.and.&
       &z(i)>ell2+ell1-dz2/3.and.z(i)<ell2+ell1+dz2/3) then
      print'(i6,a)',i
   endif
end do
!---載荷条件の節点群
print'("*NSET,NSET=saika")'
nsaika=0
do i=1,nset1
   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-dx1/3.and.x(i)<=itix+saikax/2+dx1/3) then
      nsaika=nsaika+1
      print'(i6,a)',i
   endif
end do
!#--支間中央の節点群
!print'("*NSET,NSET=tyuuou")'
!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)>=-dx/3.and.x(i)<=dx/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-dx/3.and.x(i)<=b1+dx/3) then
!      print'(i6,a)',i
!   endif
!   !---高さが中央の場合
!!   if (y(i)>h/2-dy/2.and.y(i)<h/2+dy/2.and.&
!!       &z(i)>=ell/2-saikaz/2-dz/2.and.z(i)<=ell/2+saikaz/2+dz/2.and.&
!!       &x(i)>=-dx/3.and.x(i)<=dx/3) then
!!      print'(i6,a)',i
!!   else if (y(i)>h/2-dy/2.and.y(i)<h/2+dy/2.and.&
!!       &z(i)>=ell/2-saikaz/2-dz/2.and.z(i)<=ell/2+saikaz/2+dz/2.and.&
!!       &x(i)>=b1-dx/3.and.x(i)<=b1+dx/3) then
!!      print'(i6,a)',i
!!   endif
!end do
!---境界条件
print'("*BOUNDARY")'
print'("kotei1,1,3")'
print'("kotei2,1,2")'
!---材料設定
print'("*MATERIAL,NAME=EL1")'!CLT
!print'("*ELASTIC")'!等方性
!print'(f12.4,a,f12.4)',ezz,',',poizx
print'("*ELASTIC,TYPE=ORTHO")'
call ortho(e66,exx1,eyy1,ezz1,poixy,poixz,poiyz,gyz1,&
& gxz1,gxy1,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=EL2")'!集成材
print'("*ELASTIC,TYPE=ORTHO")'
call ortho(e66,exx2,eyy2,ezz2,poixy,poixz,poiyz,gyz2,&
& gxz2,gxy2,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'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E2,MATERIAL=EL1,ORIENTATION=OR2")'
print'("*SOLID SECTION,ELSET=E3,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=E4,MATERIAL=EL1,ORIENTATION=OR2")'
print'("*SOLID SECTION,ELSET=E5,MATERIAL=EL1,ORIENTATION=OR1")'
print'("*SOLID SECTION,ELSET=Ehari,MATERIAL=EL2,ORIENTATION=OR1")'!集成材
!---載荷条件
print'("*STEP")'
print'("*STATIC")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/nsaika
!print'("*FREQUENCY")'!固有値解析
!print'("10")'!振動モード
!#--datファイルの設定
!print'("*NODE PRINT,NSET=tyuuou")'
!print'("U")'
!print'("*EL PRINT,ELSET=")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
print'("S,E")'
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
