!  木高研の試験体(半断面)を、トラス部分をB32、木床版をC3D8で作ってみる。
!  ※このプログラムで書いたinpファイルでは節点番号が共有されないので、
!    inpファイルの木部材要素の節点番号を共有させるトラス部材の
!    節点番号に書き換える必要がある。
!  ※そのため、木部材の高さ(y軸)方向の分割数は偶数にする。
         implicit real*8(a-h,o-z)
         dimension x(1000000),y(1000000),z(1000000)
!  節点座標：x(節点番号),y(節点番号),z(節点番号)
        print'("**")'
        print'("** Structure: Pratt truss")'
        print'("** Test objective: B32 and C3D8 elements.")'
        print'("** 木高研試験体(半断面)")'
        print'("**")'
! <節点番号ごとの節点座標>
! ※プラットトラスについて
        print'("** トラス部材の節点番号と座標")'
        print'("*NODE, NSET=Truss")'
!  主構高(y方向)：h(m)
!  支間長(z方向)：ell(m)
         h=380.d-3
         ell=4900.d-3
!  部材の断面(長方形へ換算した断面)
!  (1) 上弦材(x方向:bj、y方向:hj)
         bj=6.d-3
         hj=163.d-3
!  (2) 下弦材(x方向:bk、y方向:hk)
         bk=6.d-3
         hk=163.d-3
!  (3) 垂直材(x方向:bs、y方向:hs)
         bs=12.d-3
         hs=68.d-3
!  (3) 左の斜材(x方向:bsL、y方向:hsL)
         bsL=12.d-3
         hsL=68.d-3
!  (4) 右の斜材(x方向:bsR、y方向:hsR)
         bsR=12.d-3
         hsR=68.d-3
!  y方向の節点数：ny
         ny=3 ! 固定(今のところ)
!  上弦材(下弦材)の数：nk
         nk=14 ! 固定
!  上弦材(下弦材)の1本あたりの要素分割数：nz
         nz=2
!  上弦材(下弦材)1本あたりの長さ：zn
         zn=ell/nk
!
!  節点番号：nset
           nset=0
        do k=0,nk*nz*2
        do j=0,ny-1
           nset=nset+1
           x(nset)=0
           y(nset)=(h*j)/2
           z(nset)=(zn/(2*nz))*k
        print'(i6,a,f9.4,a,f9.4,a,f9.4)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
        end do
        end do
!
! ※木部材について(上部:ue,下部:sita)
        print'("** 木部材の節点番号と座標")'
        print'("*NODE, NSET=Wood_ue")'
!  床版の幅(x方向,半断面で)：bmu(m)
!  板厚(y方向)：hmu(m)
!  支間長(z方向)：ell(m)
         bmu=420.d-3
         hmu=120.d-3
!         ell=4900.d-3
!  x,y,z方向の要素分割数(偶数で)：mx,my,mz
         mx=4
         my=4
         mz=56
! 1要素のx,y,z方向の長さ：xn,yn,zn
         xm=bmu/mx
         ym=hmu/my
         zm=ell/mz
!  節点番号：nset
            nset=(nk*nz*2+1)*ny
        do k=0,mz
        do j=0,my
        do i=0,mx
           nset=nset+1
           x(nset)=xm*i
           y(nset)=ym*j-(hmu/2)
           z(nset)=zm*k
        print'(i6,a,f9.4,a,f9.4,a,f9.4)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
        end do
        end do
        end do
!
        print'("*NODE, NSET=Wood_sita")'
!  節点番号：nset
            nset=(nk*nz*2+1)*ny+((mx+1)*(my+1)*(mz+1))
        do k=0,mz
        do j=0,my
        do i=0,mx
           nset=nset+1
           x(nset)=xm*i
           y(nset)=ym*j-(hmu/2)+h
           z(nset)=zm*k
        print'(i6,a,f9.4,a,f9.4,a,f9.4)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
        end do
        end do
        end do
!
!
! <要素番号とその要素を構成する節点番号>
! ※プラットトラスについて(B32要素1分割)
!  わかりやすいように各部材に分別して表示させる
        print'("** トラス部材の要素番号と構成する節点番号")'
! (1) 上弦材
        print'("*ELEMENT, TYPE=B32, ELSET=jougen")'
          nyou=0
        do k=1,nk*nz
          nyou=nyou+1
         n1=6*k-5
         n2=6*k-2
         n3=6*k+1
        print'(i6,a,i6,a,i6,a,i6)',
     &   nyou,',',n1,',',n2,',',n3
        end do
! (2) 下弦材
        print'("*ELEMENT, TYPE=B32, ELSET=kagen")'
          nyou=nk*nz
        do k=1,nk*nz
          nyou=nyou+1
         n1=6*k-3
         n2=6*k
         n3=6*k+3
        print'(i6,a,i6,a,i6,a,i6)',
     &   nyou,',',n1,',',n2,',',n3
        end do
! (3) 垂直材
        print'("*ELEMENT, TYPE=B32, ELSET=suityoku")'
          nyou=nk*nz+nk*nz
        do k=1,nk+1
          nyou=nyou+1
         n1=6*nz*k-(6*nz-1)
         n2=6*nz*k-(6*nz-2)
         n3=6*nz*k-(6*nz-3)
        print'(i6,a,i6,a,i6,a,i6)',
     &   nyou,',',n1,',',n2,',',n3
        end do
! (4) 左の斜材
        print'("*ELEMENT, TYPE=B32, ELSET=syazaiL")'
          nyou=nk*nz+nk*nz+nk+1
        do k=1,nk/2
          nyou=nyou+1
         n1=6*nz*k-(6*nz-1)
         n2=6*nz*k-(3*nz-2)
         n3=6*nz*k+3
        print'(i6,a,i6,a,i6,a,i6)',
     &   nyou,',',n1,',',n2,',',n3
        end do
! (5) 右の斜材
        print'("*ELEMENT, TYPE=B32, ELSET=syazaiR")'
          nyou=nk*nz+nk*nz+nk+1+nk/2
        do k=nk/2+1,nk
          nyou=nyou+1
         n1=6*nz*k-(6*nz-3)
         n2=6*nz*k-(3*nz-2)
         n3=6*nz*k+1
        print'(i6,a,i6,a,i6,a,i6)',
     &   nyou,',',n1,',',n2,',',n3
        end do
!
! ※木部材について(C3D8要素)
        print'("** 木部材の要素番号と構成する節点番号")'
        print'("** トラス部材と共有させる節点番号を書き換える必要あり")'
        print'("*ELEMENT, TYPE=C3D8, ELSET=Banmen_ue")'
          nyou=nk*nz+nk*nz+nk+1+nk
        do k=1,mz
        do j=1,my
        do i=1,mx
          nyou=nyou+1
         n1=(mx+1)*(my+1)*(k-1)+(mx+1)*(j-1)+i+3*(nk*nz*2+1)
         n2=n1+1
         n3=n1+(mx+1)
         n4=n3+1
         n5=n1+(mx+1)*(my+1)
         n6=n2+(mx+1)*(my+1)
         n7=n3+(mx+1)*(my+1)
         n8=n4+(mx+1)*(my+1)
        print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',
     &    nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
        end do
        end do
        end do
!
        print'("*ELEMENT, TYPE=C3D8, ELSET=Banmen_sita")'
          nyou=nk*nz+nk*nz+nk*nz+1+nk+(mx*my*mz)
        do k=1,mz
        do j=1,my
        do i=1,mx
          nyou=nyou+1
         n1=(mx+1)*(my+1)*(k-1)+(mx+1)*(j-1)+i+3*(nk*nz*2+1)
     &      +((mx+1)*(my+1)*(mz+1))
         n2=n1+1
         n3=n1+(mx+1)
         n4=n3+1
         n5=n1+(mx+1)*(my+1)
         n6=n2+(mx+1)*(my+1)
         n7=n3+(mx+1)*(my+1)
         n8=n4+(mx+1)*(my+1)
        print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',
     &    nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
        end do
        end do
        end do
!
! <境界条件>
!  拘束節点の入力は、
!  (節点番号,拘束する最初の自由度,拘束する最後の自由度)
!  の順。
! ※トラスについて
        print'("** トラス部材の境界条件")'
! ヒンジ支承
        print'("*NSET,NSET=HingeT")'
!  拘束節点の節点番号
        print'(i6,a)',ny,','
! ローラー支承
        print'("*NSET,NSET=RollerT")'
        print'(i6,a)',(nk*nz*2+1)*ny,','
!        
! ※木部材について(下の部材のみ)
        print'("** 木部材の境界条件")'
! ヒンジ支承
        print'("*NSET,NSET=HingeW")'
!  拘束節点の節点番号
          HW=(nk*nz*2+1)*ny+((mx+1)*(my+1)*(mz+1))
     &       +(mx+1)*my+1
        do i=HW,HW+mx
        print'(i6,a)',i,','
        end do
! ローラー支承
        print'("*NSET,NSET=RollerW")'
          RW=(nk*nz*2+1)*ny+((mx+1)*(my+1)*(mz+1))
     &       +(mx+1)*(my+1)*mz+(mx+1)*my+1
        do i=RW,RW+mx
        print'(i6,a)',i,','
        end do
! 載荷部のz方向変位
!        print'("*NSET,NSET=SaikaZ")'
!        print*,nz*ny+1,","
!        do i=1,mx+1
!        print*,(nz*2+1)*ny+(mx+1)*(my+1)*mz/2+i,","
!        end do
! 木橋断面の中央切断面
        print'("*NSET,NSET=Danmen_ue")'
        do k=0,mz
        do j=1,my+1
        print'(i6,a)',(nk*nz*2+1)*ny+(mx+1)*j+(mx+1)*(my+1)*k,','
        end do
        end do
!
        print'("*NSET,NSET=Danmen_sita")'
        do k=0,mz
        do j=1,my+1
        print'(i6,a)',(nk*nz*2+1)*ny+(mx+1)*(my+1)*(mz+1)+(mx+1)*j
     &                  +(mx+1)*(my+1)*k,','
        end do
        end do
!
!  上記(NSET)の拘束節点の拘束変位
        print'("** それぞれの拘束変位")'
        print'("*BOUNDARY")'
        print'(" HingeT,      1,2")'
        print'(" HingeT,      5,6")'
        print'(" RollerT,     1,2")'
        print'(" RollerT,     5,6")'
        print'(" HingeW,      1,2")'
        print'(" RollerW,     1,2")'
!        print'(" SaikaZ,      3")'
        print'(" Danmen_ue,   1")'
        print'(" Danmen_sita, 1")'
!
!
! <材料の定義>
!  3行目にヤング率(MPa),ポアソン比
! ※トラスについて
        print'("** トラス部材に用いる材料")'
        print'("*MATERIAL,NAME=METAL")'
        print'("*ELASTIC")'
        print'(" 206.d3, 0.3")'
! ※木部材について(C3D8要素)
        print'("** 木部材に用いる材料")'
        print'("*MATERIAL,NAME=WOOD")'
! *ELASTICは等方性の意味
        print'("*ELASTIC")'
        print'(" 7.11d3, 0.4")'
!
!
! <材料と要素の関連づけ>
! ※プラットトラスについて
! (1) 上弦材
       print'("*BEAM SECTION,ELSET=jougen,MATERIAL=METAL,SECTION=RECT")'
!  断面の第1軸方向の長さ,第2軸方向の長さ
!  断面の第1軸を定義するための単位ベクトル
!        print'(i6,a)',i,','
        print*,bj,',',hj
        print*,1.,',',0.,',',0.
! (2) 下弦材
        print'("*BEAM SECTION,ELSET=kagen,MATERIAL=METAL,SECTION=RECT")'
        print*,bk,',',hk
        print*,1.,',',0.,',',0.
! (3) 垂直材
       print'("*BEAM SECTION,ELSET=suityoku,MATERIAL=METAL,SECTION=RECT"
     &  )'
        print*,bs,',',hs
        print*,1.,',',0.,',',0.
! (4) 左の斜材
      print'("*BEAM SECTION,ELSET=syazaiL,MATERIAL=METAL,SECTION=RECT")'
        print*,bsL,',',hsL
        print*,1.,',',0.,',',0.
! (5) 右の斜材
      print'("*BEAM SECTION,ELSET=syazaiR,MATERIAL=METAL,SECTION=RECT")'
        print*,bsR,',',hsR
        print*,1.,',',0.,',',0.
! ※木部材について
        print'("*SOLID SECTION,ELSET=Banmen_ue,MATERIAL=WOOD")'
        print'("*SOLID SECTION,ELSET=Banmen_sita,MATERIAL=WOOD")'
!
! <載荷条件>
! 節点ごとに異なる載荷条件を与える場合をふまえて、
! 載荷条件の種類に対して適当な名前(LASTとか)を付ける
!
! ※プラットトラスについて(中央載荷)
        print'("** トラス部材の載荷節点")'
        print'("*NSET,NSET=SaikaT")'
        print*,nk*nz*ny+1,","
! ※木部材について(中央載荷)
        print'("** 木部材の載荷節点")'
        print'("*NSET,NSET=SaikaWhasi")'
        print*,(nk*nz*2+1)*ny+(mx+1)*(my+1)*mz/2+1,","
        print*,(nk*nz*2+1)*ny+(mx+1)*(my+1)*mz/2+mx+1,","
!
        print'("*NSET,NSET=SaikaWnaka")'
        do i=2,mx
        print*,(nk*nz*2+1)*ny+(mx+1)*(my+1)*mz/2+i,","
        end do
!
!
        print'("*STEP")'
        print'("*STATIC,SOLVER=SPOOLES")'
        print'("*CLOAD")'
! 載荷節点のy(つまり2)方向に荷重を分配して載荷する。
! P/(mx+2)載荷
        write(*,*)'SaikaT,    2,',27.5d-3/(mx+2)
! P/(2*(mx+2))載荷
        write(*,*)'SaikaWhasi,2,',27.5d-3/((mx+2)*2)
! P/(4*(mx+2))載荷
        write(*,*)'SaikaWnaka,2,',27.5d-3/((mx+2)*4)
!
!        print'("**BUCKLE")'
!        print'("**10,0.1")'
        print'("*NODE PRINT,NSET=Truss")'
        print'("U")'
        print'("*NODE PRINT,NSET=Wood_ue")'
        print'("U")'
        print'("*NODE PRINT,NSET=Wood_sita")'
        print'("U")'
        print'("*EL PRINT,ELSET=jougen")'
        print'("S")'
        print'("*EL PRINT,ELSET=kagen")'
        print'("S")'
        print'("*EL PRINT,ELSET=suityoku")'
        print'("S")'
        print'("*EL PRINT,ELSET=syazaiL")'
        print'("S")'
        print'("*EL PRINT,ELSET=syazaiR")'
        print'("S")'
        print'("*EL PRINT,ELSET=Banmen_ue")'
        print'("S")'
        print'("*EL PRINT,ELSET=Banmen_sita")'
        print'("S")'
        print'("*NODE FILE,output=3D")'
        print'("U")'
        print'("*EL FILE,POSITION=AVERAGED AT NODES")'
        print'("S")'
        print'("*END STEP")'
!
!
!
      end
