! beam8p.inpの要領で2軸対称I型梁の入力データ生成のプログラム。
! 07/11/12 調整中
         implicit real*8(a-h,o-z)
         dimension x(100000),y(100000),z(100000)
! 節点座標：x(節点番号),y(節点番号),z(節点番号)
         print'("**")'
         print'("** beam8p.inpと同じデータを要素番号と節点番号の")'
         print'("** 並び順(と要素数)を変えたデータ")'
         print'("**")'
         print'("*NODE, NSET=Nall")'
        b=10.
        b=300.d-2
        b=304.8d-3 ! フランジ幅(m)
!        
        tf=20.d-2
        tf=12.7d-3 ! フランジ厚(m)
!        
        tw=2.
        tw=12.d-2
        tw=9.5d-3  ! ウェブ厚(m)
!        
        hw=10.
        hw=548.d-2
        hw=610.d-3 ! ウェブ高(m)
!        
        ell=40.
        ell=50.d0
        ell=10.d0 ! 梁の軸長(z方向)：ell(m)
! 要素分割
        nz=100 ! 軸(z)方向の要素分割
        nfx=6 ! フランジ幅(x方向)の要素分割(偶数)
        nfy=4 ! フランジ厚(y方向)の要素分割
        nwx=4 ! ウェブ厚(x方向)の要素分割(偶数)
        nwy=20 ! ウェブ高(y方向)の要素分割
        nsode=(nfx-nwx)/2 ! フランジ片袖部分の幅方向要素分割
! 静的載荷問題の場合の荷重 pload
         pload=1.d0
! 静的載荷で変位を求める問題か、座屈解析か
!        nzakutu=0 ! 片持ち梁のたわみ
!        nzakutu=1 ! 圧縮柱のオイラー座屈
         nzakutu=2 ! 梁の横ねじれ座屈
!
! そり拘束するかどうか(上下フランジの拘束の有無)
        nsori=0 !そり拘束する
!       nsori=1 !そり拘束しない
!
!
! 3種類の要素(フランジ袖、フランジ中央、ウェブ)のx,y,z方向の長さ
         zn=ell/nz ! z方向の長さ(3種類の要素共通)
         xnf=(b-tw)/(nfx-nwx) ! フランジ袖部分要素のx方向の長さ
         ynf=tf/nfy ! フランジ部分要素のy方向の長さ
         xnw=tw/nwx ! ウェブ部分要素のx方向の長さ
         ynw=hw/nwy ! ウェブ部分要素のy方向の長さ
!
!         print*,'要素数：',nx*ny*nz
!         print*,'節点数：',(nx+1)*(ny+1)*(nz+1)
! 節点番号ごとの節点座標
! 節点番号：nset
! 上フランジ
           nset=0
        do k=0,nz
        do j=0,nfy
        do i=0,nfx
           nset=nset+1
           if(i<nsode) then ! 左袖部分
            x(nset)=xnf*i
           else if(i<=nsode+nwx) then ! ウェブ直上
            x(nset)=(b-tw)/2.d0
     &              +xnw*(i-nsode)
           else ! 右袖部分
            x(nset)=(b-tw)/2.d0+tw
     &              +xnf*(i-nsode-nwx)
           end if
           y(nset)=ynf*j
           z(nset)=zn*k
        print'(i6,a,f9.3,a,f9.3,a,f9.3)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
         end do
!        print*,'nset=',nset
!
! ウェブ
        do k=0,nz
        do j=1,nwy-1 ! 上下フランジに接しない節点
        do i=0,nwx
           nset=nset+1
           x(nset)=(b-tw)/2.d0+xnw*i
           y(nset)=tf+ynw*j
           z(nset)=zn*k
        print'(i6,a,f9.3,a,f9.3,a,f9.3)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
         end do
!        print*,'nset=',nset
!
! 下フランジ
        do k=0,nz
        do j=0,nfy
        do i=0,nfx
           nset=nset+1
           if(i<nsode) then ! 左袖部分
            x(nset)=xnf*i
           else if(i<=nsode+nwx) then ! ウェブ直上
            x(nset)=(b-tw)/2.d0
     &              +xnw*(i-nsode)
           else ! 右袖部分
            x(nset)=(b-tw)/2.d0+tw
     &              +xnf*(i-nsode-nwx)
           end if
            y(nset)=tf+hw+ynf*j
            z(nset)=zn*k
        print'(i6,a,f9.3,a,f9.3,a,f9.3)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
         end do
!        print*,'nset=',nset
!
!
!
! 要素番号とその要素の8節点の節点番号
! 上フランジ
         nyou=0
        do k=1,nz
        do j=1,nfy
        do i=1,nfx
         nyou=nyou+1
         n1=(nfx+1)*(nfy+1)*(k-1)+(nfx+1)*(j-1)+i
         n2=n1+1
         n3=n1+(nfx+1)
         n4=n3+1
         n5=n1+(nfx+1)*(nfy+1)
         n6=n2+(nfx+1)*(nfy+1)
         n7=n3+(nfx+1)*(nfy+1)
         n8=n4+(nfx+1)*(nfy+1)
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         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
!
! ウェブ
        nsetf=(nfx+1)*(nfy+1)*(nz+1) !上フランジの全節点数
        nsetw=nsetf+(nwx+1)*(nwy-1)*(nz+1) !ウェブ下面の1層上までの節点数
        do k=1,nz
        do j=1,nwy
        do i=1,nwx
         nyou=nyou+1
        if(j==1) then ! 上フランジと接する要素の節点:nu
         nu1=(nfx+1)*(nfy+1)*(k-1)+(nfx+1)*nfy+nsode+i
         nu2=nu1+1
         nu5=nu1+(nfx+1)*(nfy+1)
         nu6=nu5+1
         n1=nsetf+(nwx+1)*(nwy-1)*(k-1)+(nwx+1)*(j-1)+i
         n2=n1+1
         n3=n1+(nwx+1)
         n4=n3+1
         n5=n1+(nwx+1)*(nwy-1)
         n6=n2+(nwx+1)*(nwy-1)
         n7=n3+(nwx+1)*(nwy-1)
         n8=n4+(nwx+1)*(nwy-1)
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',
     & nyou,',',nu1,',',n1,',',n5,',',nu5,',',nu2,',',n2,',',n6,',',nu6
        else if(j<nwy) then ! 上下フランジに接しない要素
         n1=nsetf+(nwx+1)*(nwy-1)*(k-1)+(nwx+1)*(j-2)+i
         n2=n1+1
         n3=n1+(nwx+1)
         n4=n3+1
         n5=n1+(nwx+1)*(nwy-1)
         n6=n2+(nwx+1)*(nwy-1)
         n7=n3+(nwx+1)*(nwy-1)
         n8=n4+(nwx+1)*(nwy-1)
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         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
        else ! 下フランジに接する要素の節点:ns
         ns3=nsetw+(nfx+1)*(nfy+1)*(k-1)+nsode+i
         ns4=ns3+1
         ns7=ns3+(nfx+1)*(nfy+1)
         ns8=ns7+1
         n1=nsetf+(nwx+1)*(nwy-1)*(k-1)+(nwx+1)*(j-2)+i
         n2=n1+1
         n3=n1+(nwx+1)
         n4=n3+1
         n5=n1+(nwx+1)*(nwy-1)
         n6=n2+(nwx+1)*(nwy-1)
         n7=n3+(nwx+1)*(nwy-1)
         n8=n4+(nwx+1)*(nwy-1)
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',
     & nyou,',',n1,',',ns3,',',ns7,',',n5,',',n2,',',ns4,',',ns8,',',n6

        end if
        end do
        end do
        end do
!
! 下フランジ
        do k=1,nz
        do j=1,nfy
        do i=1,nfx
         nyou=nyou+1
         n1=nsetw+(nfx+1)*(nfy+1)*(k-1)+(nfx+1)*(j-1)+i
         n2=n1+1
         n3=n1+(nfx+1)
         n4=n3+1
         n5=n1+(nfx+1)*(nfy+1)
         n6=n2+(nfx+1)*(nfy+1)
         n7=n3+(nfx+1)*(nfy+1)
         n8=n4+(nfx+1)*(nfy+1)
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         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
!
!
!
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の種類に対して適当な名前(FIXとか)を付ける
! 後藤メモ：FIXの意味が上記でいいかどうか確認して下さい
         print'("*NSET,NSET=FIX")'
! 拘束節点の節点番号
         if(nsori==0) then !そり拘束のときは全部拘束
          do i=1,(nfx+1)*(nfy+1) !上フランジ固定端
          print'(i6,a)',i,','
          end do
          do i=nsetf+1, nsetf+(nwx+1)*(nwy-1) !ウェブ固定端
          print'(i6,a)',i,','
          end do
          do i=nsetw+1, nsetw+(nfx+1)*(nfy+1) !下フランジ固定端
          print'(i6,a)',i,','
          end do
         else  !そり自由のときは、ウェブとウェブ直上下
          do k=1,nfy+1
          do i=nfx/2+1+(k-1)*(nfx+1)-nwx/2,nfx/2+1+(k-1)*(nfx+1)+nwx/2
                                                         !上フランジ固定端
          print'(i6,a)',i,','
          end do
          end do
          do i=nsetf+1, nsetf+(nwx+1)*(nwy-1) !ウェブ固定端
          print'(i6,a)',i,','
          end do
          do k=1,nfy+1
          do i=nsetw+nfx/2+1+(k-1)*(nfx+1)-nwx/2,
     &           nsetw+ nfx/2+1+(k-1)*(nfx+1)+nwx/2 !下フランジ固定端
          print'(i6,a)',i,','
          end do
          end do
          end if 
! 上記(NSET=FIX)の拘束節点の拘束変位
         print'("*BOUNDARY")'
         print'("FIX,1,3")'
! 節点ごとに異なる拘束条件を与える場合は、
! FIX2とかとして上記を繰り返すのかな(? 調べて)。
         print'("*MATERIAL,NAME=EL")'
         print'("*ELASTIC")'
! ヤング率とポアソン比(てことは異方性のときはどうすんのかな)
         !e=210000. ! ヤング率(MPa)
          e=21000000. ! ヤング率(1.d4*Pa) 座屈出る
          e=210000000. ! ヤング率(kPa) 座屈出る
         poi=0.3   !  ポアソン比
         print*, e,',',poi
         print'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")'
!
! 載荷条件
! 節点ごとに異なる載荷条件を与える場合をふまえて、
! 載荷条件の種類に対して適当な名前(LASTとか)を付ける
! 後藤メモ：LASTの意味が上記でいいかどうか確認して下さい
!
!
         print'("*NSET,NSET=LAST")'
!
         nhasif=nsetf-(nfx+1)*(nfy+1)+1
         do i=nhasif,nsetf !自由端上フランジ
         print'(i6,a)',i,','
         end do
         nhasiw=nsetw-(nwx+1)*(nwy-1)+1
         do i=nhasiw,nsetw !自由端ウェブ
         print'(i6,a)',i,','
         end do
         nhasi=nset-(nfx+1)*(nfy+1)+1
         do i=nhasi,nset !自由端下フランジ
         print'(i6,a)',i,','
         end do
!
         nx1ny1=(nfx+1)*(nfy+1)*2+(nwx+1)*(nwy-1)
         print'("*STEP")'
         if(nzakutu==0) then
          print'("*STATIC,SOLVER=SPOOLES")'
          print'("*CLOAD")'
          print*,'LAST,2,',pload/nx1ny1
         else if(nzakutu==1) then
         print'("*BUCKLE")'
         print*,'5,0.01'
         print'("*CLOAD")'
         print*,'LAST,3,',-1.d0/nx1ny1
           yi=b*h**3/12.d0
           pi=2.d0*asin(1.d0)
           euler=e*yi*(pi/2.d0/ell)**2
           open(8,file='euler')
           write(8,*) euler
           close(8)
         else if(nzakutu==2) then
         print'("*BUCKLE")'
         print*,'5,0.01'
         print'("*CLOAD")'
         print*,'LAST,2,',1.d0/nx1ny1
         end if
! 節点ごとに異なる載荷条件を与える場合は、
! LAST2とかとして上記を繰り返すのかな(? 調べて)。
!
         print'("*NODE PRINT,NSET=Nall")'
         print'("U")'
         print'("*EL PRINT,ELSET=Eall")'
         print'("S")'
         print'("*NODE FILE")'
         print'("U")'
         print'("*EL FILE,POSITION=AVERAGED AT NODES")'
         print'("S")'
         print'("*END STEP")'
!
!
!
         end
