!CCXでオンサイト木橋(孔無し非一体化モデル)のリラクセーションの1/4解析
! プレストレス-->載荷-->プレストレスだけ減少 でとずれがどうなるのかを計算。
! 現場で計測している緊張力がどこまで減っていても安全といえるかを見積る。
implicit none
!#---入力
real    :: b1=840.d-3/2       ! m   :木材全体の幅(x)/2
real    :: h1=120.d-3         ! m   :一本の木材の高さ=幅(y or z)
real    :: h2=260.d-3         ! m   :上下木材の間隔(y)
real    :: b2=9.0d-3          ! m   :鋼板の厚さ(x)
real    :: ell=7300.d-3/2     ! m   :長さ(z)/2
real    :: koteiy=150.d-3     ! m   :端から境界条件までの距離(z)
real    :: py=70.d-3/4        ! MN  :荷重/4(二等林道橋の後輪荷重55kN*(1+木橋の衝撃荷重0.25))
real    :: px0=70.d-3         ! MN  :導入時のプレストレス力
real    :: px=40.d-3          ! MN  :緊張力低下後のプレストレス力
real    :: saikabanx=500d-3/2 ! m   :載荷板の長さ(x)
real    :: saikabanz=200d-3/2 ! m   :載荷板の横幅(z)
real    :: henni=170d-3        ! m   :中央から変位計までの距離(z)
real    :: anka=100.d-3       ! m   :アンカープレートの一辺の長さ
real    :: hazip=350.d-3      ! m   :端からPC鋼棒までの距離(z)
real    :: ellp=700.d-3       ! m   :PC鋼棒の間隔
integer :: np=10               !     :PC鋼棒の本数(上下合わせて)
integer :: nstep1=1           !     :プレストレスの出力ステップ数
integer :: nstep2=1           !     :荷重の出力ステップ数
integer :: nstep3=30          !     :緊張力低下時の出力ステップ数
integer :: nx1=8              !偶数 :一本の木材の幅員方向(x)分割数
integer :: nx2=4              !     :鋼板の厚さ方向(x)分割数
integer :: ny1=8              !     :上下木材の高さ方向(y)分割数
integer :: ny2=16             !     :上下木材間隔の高さ方向(y)分割数
integer :: nz=128             !     :橋軸方向(z)分割数
real    :: ezz=7.694d3!6.921d3        ! MPa :木材のz軸方向ヤング率
real    :: exx=7.694d3/25.0d0 ! MPa :木材のx軸方向ヤング率
real    :: eyy=7.694d3/25.0d0 ! MPa :木材のy軸方向ヤング率
real    :: poizx=0.4d0        !     :木材のポアソン比
real    :: ekou=206.0d3       ! MN  :鋼材のヤング率
real    :: poikou=0.3         !     :鋼材のポアソン比
real    :: mukou=0.4           !     :鋼材と木材の摩擦係数
real    :: mumoku=0.4          !     :木材同士の摩擦係数
!---
real    :: gxy,gxz,gyz,poizy,poixy,poixz,poiyx,poiyz
real    :: x(1000000),y(1000000),z(1000000),e66(6,6),b(3,3)
integer :: mesh1(300000,8),mesh2(300000,8)
real    :: dx1,dx2,dy1,dy2,dz,ts,sakaix
integer :: i,j,k,l,nset,nset1,nset2,nset3,nyou,nyou1,nyou2,nyou3
integer :: nsaika1,nsaika2,npre
!sakaix=b2+b1*2.d0/7! m :鋼板から境界までの距離(x)
!---せん断弾性係数
gxy=ezz/15
gxz=gxy
gyz=gxy
!---ポアソン比
poizy=poizx
poixy=poizx*exx/ezz
poixz=poixy
poiyx=poixy
poiyz=poixy
!---要素長さ
dx1=h1/nx1
dx2=b2/nx2
dy1=h1/ny1
dy2=h2/ny2
dz=ell/nz
!---鋼板の節点座標
nset=0
do k=0,nz
   !---鋼板上部
   do j=0,ny1
      do i=0,nx2
         nset=nset+1
         x(nset)=dx2*i
         y(nset)=dy1*j
         z(nset)=dz*k
      enddo
   enddo
   !---鋼板中部
   do j=1,ny2
      do i=0,nx2
         nset=nset+1
         x(nset)=dx2*i
         y(nset)=h1+dy2*j
         z(nset)=dz*k
      enddo
   enddo
   !---鋼板下部
   do j=1,ny1
      do i=0,nx2
         nset=nset+1
         x(nset)=dx2*i
         y(nset)=h1+h2+dy1*j
         z(nset)=dz*k
      enddo
   enddo
enddo
nset1=nset!鋼材の総接点数
!---木材の節点座標
do l=0,int(b1/h1)-1
   !---木材上部
   do k=0,nz
      do j=0,ny1
         do i=0,nx1
            nset=nset+1
            x(nset)=b2+h1*l+dx1*i
            y(nset)=dy1*j
            z(nset)=dz*k
         enddo
      enddo
   enddo
   !---木材下部
   do k=0,nz
      do j=0,ny1
         do i=0,nx1
            nset=nset+1
            x(nset)=b2+h1*l+dx1*i
            y(nset)=h1+h2+dy1*j
            z(nset)=dz*k
         enddo
      enddo
   enddo
enddo
nset2=(nset-nset1)/int(b1/h1)/2!一本の木材の総接点数
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   l=int(b1/h1)
   !---木材上部
   do k=0,nz
      do j=0,ny1
         do i=0,int(nx1/2)
            nset=nset+1
            x(nset)=b2+h1*l+dx1*i
            y(nset)=dy1*j
            z(nset)=dz*k
         enddo
      enddo
   enddo
   !---木材下部
   do k=0,nz
      do j=0,ny1
         do i=0,int(nx1/2)
            nset=nset+1
            x(nset)=b2+h1*l+dx1*i
            y(nset)=h1+h2+dy1*j
            z(nset)=dz*k
         enddo
      enddo
   enddo
   nset3=nset/2!半分の木材の総接点数
endif
!---鋼板の要素を構成する節点
nyou=0
do k=0,nz-1
   do j=0,(ny1*2+ny2)-1
      do i=1,nx2
         nyou=nyou+1
         mesh1(nyou,1)=(nx2+1)*(ny1*2+ny2+1)*k+(nx2+1)*j+i
         mesh1(nyou,2)=mesh1(nyou,1)+(nx2+1)
         mesh1(nyou,3)=mesh1(nyou,2)+(nx2+1)*(ny1*2+ny2+1)
         mesh1(nyou,4)=mesh1(nyou,1)+(nx2+1)*(ny1*2+ny2+1)
         mesh1(nyou,5)=mesh1(nyou,1)+1
         mesh1(nyou,6)=mesh1(nyou,5)+(nx2+1)
         mesh1(nyou,7)=mesh1(nyou,6)+(nx2+1)*(ny1*2+ny2+1)
         mesh1(nyou,8)=mesh1(nyou,5)+(nx2+1)*(ny1*2+ny2+1)
      enddo
   enddo
enddo
nyou1=nyou!鋼板の総要素数
!---木材の要素を構成する節点
nyou=0
do l=0,int(b1/h1)-1
   !---木材上部
   do k=0,nz-1
      do j=0,ny1-1
         do i=1,nx1
            nyou=nyou+1
            mesh2(nyou,1)=nset1+nset2*2*l+(nx1+1)*(ny1+1)*k+(nx1+1)*j+i
            mesh2(nyou,2)=mesh2(nyou,1)+nx1+1
            mesh2(nyou,3)=mesh2(nyou,2)+(nx1+1)*(ny1+1)
            mesh2(nyou,4)=mesh2(nyou,1)+(nx1+1)*(ny1+1)
            mesh2(nyou,5)=mesh2(nyou,1)+1
            mesh2(nyou,6)=mesh2(nyou,5)+nx1+1
            mesh2(nyou,7)=mesh2(nyou,6)+(nx1+1)*(ny1+1)
            mesh2(nyou,8)=mesh2(nyou,5)+(nx1+1)*(ny1+1)
         enddo
      enddo
   enddo
   !---木材下部
   do k=0,nz-1
      do j=0,ny1-1
         do i=1,nx1
            nyou=nyou+1
            mesh2(nyou,1)=nset1+nset2+nset2*2*l+(nx1+1)*(ny1+1)*k+(nx1+1)*j+i
            mesh2(nyou,2)=mesh2(nyou,1)+nx1+1
            mesh2(nyou,3)=mesh2(nyou,2)+(nx1+1)*(ny1+1)
            mesh2(nyou,4)=mesh2(nyou,1)+(nx1+1)*(ny1+1)
            mesh2(nyou,5)=mesh2(nyou,1)+1
            mesh2(nyou,6)=mesh2(nyou,5)+nx1+1
            mesh2(nyou,7)=mesh2(nyou,6)+(nx1+1)*(ny1+1)
            mesh2(nyou,8)=mesh2(nyou,5)+(nx1+1)*(ny1+1)
         enddo
      enddo
   enddo
enddo
nyou2=nyou/int(b1/h1)/2!一本の木材の総要素数
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   l=int(b1/h1)
   !---木材上部
   do k=0,nz-1
      do j=0,ny1-1
         do i=1,int(nx1/2)
            nyou=nyou+1
            mesh2(nyou,1)=nset1+nset2*2*l+(nx1/2+1)*(ny1+1)*k+(nx1/2+1)*j+i
            mesh2(nyou,2)=mesh2(nyou,1)+nx1/2+1
            mesh2(nyou,3)=mesh2(nyou,2)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,4)=mesh2(nyou,1)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,5)=mesh2(nyou,1)+1
            mesh2(nyou,6)=mesh2(nyou,5)+nx1/2+1
            mesh2(nyou,7)=mesh2(nyou,6)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,8)=mesh2(nyou,5)+(nx1/2+1)*(ny1+1)
         enddo
      enddo
   enddo
   !---木材下部
   do k=0,nz-1
      do j=0,ny1-1
         do i=1,int(nx1/2)
            nyou=nyou+1
            mesh2(nyou,1)=nset1+nset2*2*l+(nx1/2+1)*(ny1+1)*((nz+1)+k)+(nx1/2+1)*j+i
            mesh2(nyou,2)=mesh2(nyou,1)+nx1/2+1
            mesh2(nyou,3)=mesh2(nyou,2)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,4)=mesh2(nyou,1)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,5)=mesh2(nyou,1)+1
            mesh2(nyou,6)=mesh2(nyou,5)+nx1/2+1
            mesh2(nyou,7)=mesh2(nyou,6)+(nx1/2+1)*(ny1+1)
            mesh2(nyou,8)=mesh2(nyou,5)+(nx1/2+1)*(ny1+1)
         enddo
      enddo
   enddo
   nyou3=(nyou-nyou2*int(b1/h1)*2)/2!半分の木材の総要素数
endif
!---節点座標の出力
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)
enddo
!---鋼板の要素の出力
print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")'
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)
end do
!---木材の要素の出力
do j=0,int(b1/h1)*2-1
   print'(A,i1)',"*ELEMENT, TYPE=C3D8, ELSET=Emoku",j+1
   do i=nyou2*j+1,nyou2*(j+1)
      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
enddo
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   l=int(b1/h1)*2
   print'(A,i1)',"*ELEMENT, TYPE=C3D8, ELSET=Emoku",l+1
   do i=nyou2*l+1,nyou2*l+nyou3
      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,i1)',"*ELEMENT, TYPE=C3D8, ELSET=Emoku",l+2
   do i=nyou2*l+nyou3+1,nyou2*l+nyou3*2
      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
endif
!---鋼板のマスター面群
print'(A)',"*SURFACE,NAME=SMkou1"
do k=0,nz-1
   !---鋼板上部
   do j=0,ny1-1
      l=nx2*(ny1*2+ny2)*k+nx2*j+nx2
      print'(i6,A)',l,',S2'
   enddo
   !---鋼板中部(上1/4)
   do j=0,int(ny2/4)-1
      l=nx2*(ny1*2+ny2)*k+nx2*(ny1+j)+nx2
      print'(i6,A)',l,',S2'
   enddo
enddo
print'(A)',"*SURFACE,NAME=SMkou2"
do k=0,nz-1
   !---鋼板中部(下1/4)
   do j=ny2-int(ny2/4),ny2-1
      l=nx2*(ny1*2+ny2)*k+nx2*(ny1+j)+nx2
      print'(i6,A)',l,',S2'
   enddo
   !---鋼板下部
   do j=0,ny1-1
      l=nx2*(ny1*2+ny2)*k+nx2*(ny1+ny2+j)+nx2
      print'(i6,A)',l,',S2'
   enddo
enddo
!---木材のスレーブ面群とマスター面群
do l=0,int(b1/h1)*2-1
   print'(A,i1)',"*SURFACE,NAME=SSmoku",l+1
   do k=0,nz-1
      do j=0,ny1-1
         i=nyou1+nyou2*l+(nx1*ny1)*k+nx1*j+1
         print'(i6,A)',i,',S1'
      enddo
   enddo
   print'(A,i1)',"*SURFACE,NAME=SMmoku",l+1
   do k=0,nz-1
      do j=0,ny1-1
         i=nyou1+nyou2*l+(nx1*ny1)*k+nx1*j+nx1
         print'(i6,A)',i,',S2'
      enddo
   enddo
enddo
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   l=int(b1/h1)*2
   print'(A,i1)',"*SURFACE,NAME=SSmoku",l+1
   do k=0,nz-1
      do j=0,ny1-1
         i=nyou1+nyou2*l+(nx1/2*ny1)*k+nx1/2*j+1
         print'(i6,A)',i,',S1'
      enddo
   enddo
   print'(A,i1)',"*SURFACE,NAME=SSmoku",l+2
   do k=0,nz-1
      do j=0,ny1-1
         i=nyou1+nyou2*l+nyou3+(nx1/2*ny1)*k+nx1/2*j+1
         print'(i6,A)',i,',S1'
      enddo
   enddo
endif
!---境界条件の節点群
print'("*NSET,NSET=koteiy")'!支持部
do i=1,nset
   if (y(i)>h1*2+h2-dy1/3.and.y(i)<h1*2+h2+dy1/3.and.&
       &z(i)>koteiy-dz/2.and.z(i)<koteiy+dz/2) then
      print'(i6,a)',i,','
   endif
end do
print'("*NSET,NSET=koteix")'!x軸対称面
do i=1,nset
   if (x(i)>b1+b2-dx1/3.and.x(i)<b1+b2+dy1/3) then
      print'(i6,a)',i,','
   endif
end do
print'("*NSET,NSET=koteiz")'!z軸対称面
do i=1,nset
   if (z(i)>ell-dz/3.and.x(i)<ell+dz/3) then
      print'(i6,a)',i,','
   endif
end do
!---荷重条件の節点群
print'("*NSET,NSET=saika1")'!載荷板に載荷(角材中央の節点)
nsaika1=0
nsaika2=0
do j=int(b1/h1-1.d-3)-int(saikabanx/h1-1.d-3),int(b1/h1+1.d-3)-1
   do i=nset1+nset2*2*j+1,nset1+nset2*2*(j+1)!角材のみに載荷
      if (x(i)>b2+b1-saikabanx+dx1/3.and.x(i)<b1+b2-dx1/3.and.&
          &x(i)>b2+h1*j+dx1/3.and.x(i)<b2+h1*(j+1)-dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika1=nsaika1+1
         print'(i6,a)',i,','
      end if
   end do
enddo
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   j=int(b1/h1-1.d-3) 
   do i=nset1+nset2*2*j+1,nset1+nset2*2*(j+1)!角材のみに載荷
      if (x(i)>b2+b1-saikabanx+dx1/3.and.x(i)<b1+b2-dx1/3.and.&
          &x(i)>b2+h1*j+dx1/3.and.x(i)<b2+h1*(j+1.d0/2)-dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika1=nsaika1+1
         print'(i6,a)',i,','
      end if
   end do
endif
print'("*NSET,NSET=saika2")'!載荷板に載荷(角材端の節点)
do j=int(b1/h1-1.d-3)-int(saikabanx/h1-1.d-3),int(b1/h1+1.d-3)-1
   do i=nset1+nset2*2*j+1,nset1+nset2*2*(j+1)!角材のみに載荷
      if (x(i)>b2+b1-saikabanx-dx1/2.and.x(i)<b2+b1-saikabanx+dx1/2.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika2=nsaika2+1
         print'(i6,a)',i,','
      elseif (x(i)>b2+b1-saikabanx-dx1/3.and.x(i)<b2+b1+dx1/3.and.&
          &x(i)>b2+h1*j-dx1/3.and.x(i)<b2+h1*j+dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika2=nsaika2+1
         print'(i6,a)',i,','
      else if (x(i)>b2+b1-saikabanx-dx1/3.and.x(i)<b2+b1+dx1/3.and.&
          &x(i)>b2+h1*(j+1)-dx1/3.and.x(i)<b2+h1*(j+1)+dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika2=nsaika2+1
         print'(i6,a)',i,','
      end if
   end do
enddo
if (mod(b1,h1)>h1/2-1.d-3.and.mod(b1,h1)<h1/2+1.d-3) then!木材の本数が奇数の時
   j=int(b1/h1-1.d-3) 
   do i=nset1+nset2*2*j+1,nset1+nset2*2*(j+1)!角材のみに載荷
      if (x(i)>b2+b1-saikabanx+dx1/3.and.x(i)<b1+b2-dx1/3.and.&
          &x(i)>b2+h1*j-dx1/3.and.x(i)<b2+h1*j+dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika2=nsaika2+1
         print'(i6,a)',i,','
      else if (x(i)>b2+b1-saikabanx-dx1/3.and.x(i)<b2+b1+dx1/3.and.&
          &x(i)>b2+h1*(j+1.d0/2)-dx1/3.and.x(i)<b2+h1*(j+1.d0/2)+dx1/3.and.&
          &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
          &z(i)>ell-saikabanz-dz/3.and.z(i)<ell+dz/3) then
         nsaika2=nsaika2+1
         print'(i6,a)',i,','
      end if
   end do
endif
print'(A)',"*NSET,NSET=pre"!プレストレス
npre=0
do i=1,nset
   if (x(i)>-dx2/3.and.x(i)<dx2/3.and.&!上のプレストレス
       &y(i)>h1/2-anka/2-dy1/3.and.y(i)<h1/2+anka/2+dy1/3.and.&
       &z(i)>anka/2.and.z(i)<ell-anka/2) then
      npre=npre+1
      print'(i6,a)',i,','
   end if
   if (x(i)>-dx2/3.and.x(i)<dx2/3.and.&!下のプレストレス
       &y(i)>h1+h2+h1/2-anka/2-dy1/3.and.y(i)<h1+h2+h1/2+anka/2+dy1/3.and.&
       &z(i)>anka/2.and.z(i)<ell-anka/2) then
      npre=npre+1
      print'(i6,a)',i,','
   end if
end do
print'("*NSET,NSET=henni")'!変位計の断面
do i=1,nset1!鋼板
   if (x(i)>-dx2/3.and.x(i)<b2+dx1/3.and.&
       &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>ell-dz/3.and.z(i)<ell+dz/3) then!鋼板は梁中央に変位計
      print'(i6,a)',i,','
   end if
end do
do i=nset1,nset!木材
   if (x(i)>b2-dx2/3.and.x(i)<b1+b2+dx1/3.and.&
       &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>ell-henni-dz/3.and.z(i)<ell-henni+dz/3) then!鋼板は梁中央に変位計
      print'(i6,a)',i,','
   end if
end do
print'("*NSET,NSET=jouen")'!梁中央の断面の上縁
do i=1,nset
   if (x(i)>-dx2/3.and.x(i)<b1+b2+dx1/3.and.&
       &y(i)>-dy1/3.and.y(i)<dy1/3.and.&
       &z(i)>ell-dz/3.and.z(i)<ell+dz/3) then
      print'(i6,a)',i,','
   end if
end do
print'("*NSET,NSET=kaen")'!梁中央の断面の下縁(木材のみ)
do i=nset1+1,nset
   if (x(i)>b2-dx2/3.and.x(i)<b1+b2+dx1/3.and.&
       &y(i)>h1*2.d0+h2-dy1/3.and.y(i)<h1*2.d0+h2+dy1/3.and.&
       &z(i)>ell-dz/3.and.z(i)<ell+dz/3) then
      print'(i6,a)',i,','
   end if
end do
!---境界条件
print'("*BOUNDARY")'
print'("koteiy,2")'
print'("koteix,1")'
print'("koteiz,3")'
!#--鋼板と木材の接触および摩擦の設定
do i=1,2
   print'(A,i1)',"*CONTACT PAIR,INTERACTION=SI",i
   print'(A,i1,A,i1)',"SSmoku",i,",SMkou",i!スレーブ面,マスター面
   print'(A,i1)',"*SURFACE INTERACTION,NAME=SI",i
   print'("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=EXPONENTIAL")'
   print'("1.e-4,1.e-0")'!接触と判定する距離(m),距離0の時の応力?
   print'("*FRICTION")'
   print'(1pe13.5,a,1pe13.5)',mukou,",",ezz!摩擦係数,ペナルティー法のk/|t|(収束する範囲で大きめがいい)
enddo
!#--木材同士の接触および摩擦の設定
do i=1,int(b1/h1-1.d-3)*2
   print'(A,i1)',"*CONTACT PAIR,INTERACTION=SI",i+2
   print'(A,i1,A,i1)',"SSmoku",i+2,",SMmoku",i!スレーブ面,マスター面
   print'(A,i1)',"*SURFACE INTERACTION,NAME=SI",i+2
   print'("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=EXPONENTIAL")'
   print'("1.e-4,1.e-0")'!接触と判定する距離(m),距離0の時の応力?
   print'("*FRICTION")'
   print'(1pe13.5,a,1pe13.5)',mumoku,",",ezz!摩擦係数,ペナルティー法のk/|t|(収束する範囲で大きめがいい)
enddo
!---木材の材料設定
print'("*MATERIAL,NAME=ELmoku")'
print'("*ELASTIC,TYPE=ORTHOROPIC")'
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,a)',e66(4,4)
!---鋼材の材料設定
print'("*MATERIAL,NAME=ELkou")'
print'("*ELASTIC")'
print'(f12.4,a,f12.4)',ekou,',',poikou
do i=1,(int(b1/h1-1.d-3)+1)*2
   print'(A,i1,A)',"*SOLID SECTION,ELSET=Emoku",i,",MATERIAL=ELmoku"
enddo
print'("*SOLID SECTION,ELSET=Ekou,MATERIAL=ELkou")'
!---プレストレスの出力ステップ
print'("*TIME POINTS,NAME=T1")'
do i=1,nstep1
   ts=real(i)/nstep1
   print'(f12.4,a)',ts,','
enddo
!---荷重の出力ステップ
print'("*TIME POINTS,NAME=T2")'
do i=1,nstep2
   ts=real(i)/nstep2
   print'(f12.4,a)',ts,','
enddo
!---緊張力低下時の出力ステップ
print'("*TIME POINTS,NAME=T3")'
do i=1,nstep3
   ts=real(i)/nstep3
   print'(f12.4,a)',ts,','
enddo
!---載荷条件(プレストレス)
print'("*STEP,NLGEOM")'
print'("*STATIC")'
!print'("1.0,1.0")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','pre',',',1,',',px0*np/npre
print'("*NODE FILE,TIMEPOINTS=T1")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T1")'
print'("S")'
print'("*END STEP")'
!---載荷条件(荷重)
print'("*STEP,NLGEOM")'
print'("*STATIC")'
!print'("0.01,1.0")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika1',',',2,',',py/(nsaika1+nsaika2/2)
print'(A,a,i6,a,1pe13.5)','saika2',',',2,',',py/2/(nsaika1+nsaika2/2)
!#--出力の設定
print'("*NODE PRINT,NSET=henni,TIMEPOINTS=T2")'!変位計のたわみ
print'("U")'
print'("*NODE PRINT,NSET=jouen,TIMEPOINTS=T2")'!梁中央上部のたわみ
print'("U")'
!print'("*NODE PRINT,NSET=kaen,TIMEPOINTS=T2")'!梁中央下縁のたわみ
!print'("U")'
!print'("*NODE PRINT,NSET=koteiz,TIMEPOINTS=T2")'!梁中央の断面のたわみ
!print'("U")'
print'("*NODE FILE,TIMEPOINTS=T2")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T2")'
print'("S")'
print'("*END STEP")'
!---載荷条件(プレストレスの減少)
print'("*STEP,NLGEOM")'
print'("*STATIC")'
!print'("0.01,1.0")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','pre',',',1,',',(px-px0)*np/npre
!#--出力の設定
print'("*NODE PRINT,NSET=henni,TIMEPOINTS=T3")'!変位計のたわみ
print'("U")'
print'("*NODE PRINT,NSET=jouen,TIMEPOINTS=T3")'!梁中央上部のたわみ
print'("U")'
!print'("*NODE PRINT,NSET=kaen,TIMEPOINTS=T3")'!梁中央下縁のたわみ
!print'("U")'
!print'("*NODE PRINT,NSET=koteiz,TIMEPOINTS=T3")'!梁中央の断面のたわみ
!print'("U")'
print'("*NODE FILE,TIMEPOINTS=T3")'
print'("U")'
print'("*EL FILE,TIMEPOINTS=T3")'
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)
!dimension 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



