! L型綱(7型綱)の節点、エレメントのinpファイルをだしてくれるプログラム
! 07/11/28
! 節点座標：x(節点番号),y(節点番号),z(節点番号)
         implicit real*8(a-h,o-z)
         dimension x(10000000),y(10000000),z(1000000)
         print'("**")'
         print'("** L型綱の節点とエレメントをだしてくれる")'
         print'("** プログラム")'
         print'("**")'
         print'("*NODE, NSET=Nall")'
        b=0.007 ! フランジ幅(m)
!        
        tf=0.1  ! フランジ厚さ(m)
!        
        tw=0.006! ウェブ厚(m)
!        
        hw=0.5 ! ウェブ高(m)
!        
        ell=1 ! 梁の軸長(z方向)：ell(m)
! 要素分割
        nz=2 ! 軸(z)方向の要素分割
        nfx=4 ! フランジ幅(x方向)の要素分割(偶数)
        nfy=2 ! フランジ厚(y方向)の要素分割
        nwx=2 ! ウェブ厚(x方向)の要素分割(偶数)
        nwy=2 ! ウェブ高(y方向)の要素分割
        nsode=(nfx-nwx) ! フランジ袖部分の幅方向要素分割
! 静的載荷問題の場合の荷重 pload
         pload=1.d0
! 静的載荷で変位を求める問題か、座屈解析か
      !  nzakutu=0 !静的縮み
        nzakutu=1 !座屈        
! そり拘束するかどうか(上下フランジの拘束の有無)
        nsori=0 !そり拘束する
!       nsori=1 !そり拘束しない
!
!
! 3種類の要素(7型上部分、7型の袖以外の部分、7型の奥行き)のx,y,z方向の長さ
         zn=ell/nz ! z方向の長さ(要素共通)
         xnf=(b-tw)/(nfx-nwx) ! 7型の袖部分要素のx方向の長さ
         ynf=tf/nfy ! 7型の上の部分要素のy方向の長さ
         xnw=tw/nwx ! 7型の袖以外の部分要素のx方向の長さ
         ynw=hw/nwy ! 7型の袖以外の部分要素のy方向の長さ
!
!         print*,'要素数：',nx*ny*nz
!         print*,'節点数：',(nx+1)*(ny+1)*(nz+1)
! 節点番号ごとの節点座標
! 節点番号：nset
! 7型の上の部分(フランジ?部分)
           nset=0
        do k=0,nz
        do j=0,nfy
        do i=0,nfx
           nset=nset+1
           if(i<nsode) then ! 7型の左袖部分
            x(nset)=xnf*i
           else if(i<=nsode+nwx) then ! 右の棒(ウェブ?)の直上
            x(nset)=(b-tw)
     &              +xnw*(i-nsode)
           end if
           y(nset)=ynf*j
           z(nset)=zn*k
        print'(i6,a,1pd11.3,a,1pd11.3,a,1pd11.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 ! 上フランジに接しない節点
        do i=0,nwx
           nset=nset+1
           x(nset)=(b-tw)+xnw*i
           y(nset)=tf+ynw*j
           z(nset)=zn*k
        print'(i6,a,1pd11.3,a,1pd11.3,a,1pd11.3)',
     &         nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
         end do
!        print*,'nset=',nset
!
!
! 要素番号とその要素の8節点の節点番号
! 7型の上の部分(フランジ?)
         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) !7型上の部分の全節点数
        nsetw=nsetf+(nwx+1)*(nwy)*(nz+1) !全節点数
        do k=1,nz
        do j=1,nwy
        do i=1,nwx
         nyou=nyou+1
        if(j==1) then ! 7型上の部分と接する要素の節点: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)*(k-1)+(nwx+1)*(j-1)+i
         n2=n1+1
         n3=n1+(nwx+1) 
         n4=n3+1
         n5=n1+(nwx+1)*(nwy)
         n6=n2+(nwx+1)*(nwy)
         n7=n3+(nwx+1)+(nwy)
         n8=n4+(nwx+1)+(nwy)
         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>=2) then ! 7型上の部分に接しない要素
         n1=nsetf+(nwx+1)*(nwy)*(k-1)+(nwx+1)*(j-2)+i
         n2=n1+1
         n3=n1+(nwx+1)
         n4=n3+1
         n5=n1+(nwx+1)*(nwy)
         n6=n2+(nwx+1)*(nwy)
         n7=n3+(nwx+1)*(nwy)
         n8=n4+(nwx+1)*(nwy)
         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 if
      end do
      end do
      end do

!
!
!
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の名前をつける。
         print'("*NSET,NSET=FIX1")'
! その拘束節点の節点番号
         print'(i6,a,)',nfx+1,','!手前の面のもっとも右端の点の拘束
          !そり自由
! 上記(NSET=FIX1)の拘束節点の拘束変位
         print'("*BOUNDARY")'
         print'("FIX1,1,3")'
! 2番めの拘束節点
         print'("*NSET,NSET=FIX2")'
! その拘束節点の節点番号
!
!手前面の右端の上の点以外の全ての点を辺として拘束 !
! フランジの部分                                  !
      do i=2,(nfy+1)                              !
           print'(i6,a,)',(nfx+1)*i,','           !
        end do                                    !
! ウェブの部分                                    !
         do k=1,nwy                               !
        print'(i6,a,)',nsetf+(nwy+1)*k,','        !
        end do                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!手前面の全てを拘束!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!フランジ部分                                   
       !do i=1,(nfx+1)*(nfy+1)                 
       !print'(i6,a,)',i,','
!ウェブ部分
       !do k =1,(nwx+1)*nwy
       ! print'(i6,a,)',nsetf+k,','
!        
! 上記(NSET=FIX2)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX2,1")'
         print'("FIx2,3")'
! 拘束条件の名前を付ける
          print'("*NSET,NSET=FIX3")'
! その拘束節点の節点番号
         print'(i6,a,)',nsetf-(nfx+1)*nfy,','
! 上記(NSET=FIX3)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX3,1,2")'
! 拘束節点の名前を付ける
         print'("*NSET,NSET=FIX4")'
! 拘束節点の節点番号
!奥面の橋の一辺を拘束(1番上の点以外)
! フランジ部分
        do i=1,(nfy)
         print'(i6,a,)',nsetf-(nfx+1)*nfy+(nfx+1)*i,','
        end do
!ウェブの部分
        do k=0,nwy-1
         print'(i6,a,)',nsetf+nwy*(nwx+1)*(nz+1)-(nwx+1)*k,','
         end do
!奥面の全てを拘束
         !フランジの部分
        ! do i=i,(nfx+1)*(nfy+1)
        ! print'(i6,a,)',(nfx+1)*(nfy+1)*(nz)+i,','
         !ウェブの部分
         !do k=0,(nwx+1)*nz
        ! print'(i6,a,)',nsetf+(nwx+1)*(nwy+1)*nz+i,','
! 上記(NSET=FIX4)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX4,1")'
        print'("*MATERIAL,NAME=EL")'
         print'("*ELASTIC")'
! ヤング率とポアソン比(てことは異方性のときはどうすんのかな)
         !e=210000. ! ヤング率(MPa)
          e=203000. ! ヤング率(1.d4*Pa) 座屈出る
          e=e*10E3 ! ヤング率(kPa) 座屈出る
         poi=0.3   !  ポアソン比
         print*, e,',',poi
         print'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")'
!
! 載荷条件
! 節点ごとに異なる載荷条件を与える場合をふまえて、
! 載荷条件の種類に対して適当な名前(LASTとか)を付ける
! 後藤メモ：LASTの意味が上記でいいかどうか確認して下さい
!
!
         print'("*NSET,NSET=LAST")'
       ! フランジ上載荷の節点番号     
         do i=0,nfy
         print'(i6,a,)',nsetf-(nfx+1)*nfy+(nfx+1)*i,','
         end do
       ! ウェブ上の載荷の節点番号
        do k=0,nwy-1
          print'(i6,a,)',nsetf+nwy*(nwx+1)*(nz+1)-(nwx+1)*k,','
          end do
      !
      !
      !   
      !   
         print'("*STEP")'
         if(nzakutu==0) then
          print'("*STATIC,SOLVER=SPOOLES")'
          print'("*CLOAD")'
          print*,'LAST,3,',-1.d0/(nfy+nwy+1)*2
         else if(nzakutu==1) then
         print'("*BUCKLE")'
         print*,'5,0.01'
         print'("*CLOAD")'
         print*,'LAST,3,',-10.d0/(nfy+nwy+1)
           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)
         end if
! 
!
!
!
         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
