! 直方体のエレメント、節点をだしてくれるプログラム
! 09/12/10
! 節点座標：x(節点番号),y(節点番号),z(節点番号)
         implicit real*8(a-h,o-z)
         dimension x(10000000),y(10000000),z(1000000)
         print'("**")'
         print'("** 直方体に荷重をかけるプログラム")'
         print'("** プログラム")'
         print'("**")'
         print'("*NODE, NSET=Nall")'
        b=0.006 ! 直方体の横(x方向)(m)
        tf=0.06  ! 直方体の縦(y方向)(m)
        ell=1. ! 直方体の奥行き(z方向)(m)
        bairitu=1.d2
        bairitu=1.d0
        t=t*bairitu
        tf=tf*bairitu
        ell=ell*bairitu
! 要素分割
        nz=10 ! 軸(z)方向の分割
        nfx=6! 幅(x)方向の分割
        nfy=6 ! 厚さ(y)方向の分割
! 静的載荷問題の場合の荷重 pload
         pload=1.d0
! 静的載荷で変位を求める問題か、座屈解析か
      !  nzakutu=0 !静的縮み
        nzakutu=1 !座屈        
! そり拘束するかどうか(上下フランジの拘束の有無)
        ! nsori=0 !そり拘束する
        nsori=1 !そり拘束しない
!
!
! 3種類の要素の1要素分のx,y,z方向の長さ
         zn=ell/nz ! z方向の長さ
         xnf=b/nfx ! x方向の長さ
         ynf=tf/nfy ! 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
            x(nset)=xnf*i
           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
!
!
!
! 要素番号とその要素の8節点の節点番号
         print'("*ELEMENT, TYPE=C3D8, ELSET=Eall")'
         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'(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'("*NSET,NSET=FIX1")'
! その拘束節点の節点番号
         print'(i6,a,)',nfx+1,','
! 上記(NSET=FIX1)の拘束節点の拘束変位
         print'("*BOUNDARY")'
         print'("FIX1,1,3")'
! 2番めの拘束節点
         print'("*NSET,NSET=FIX2")'
! その拘束節点の節点番号         
         !
        ! do i=1, (nfx+1)*(nfy+1)
        ! print'(i6,a,)',i,','  !手前の端面をすべて拘束
        ! end do
         !
       do i=2,(nfy+1)
            print'(i6,a,)',(nfx+1)*i,','  !縁部の1辺を拘束の場合
           ! print'(i6,a,)',(nfx+1)*i-nfx/2,',' !中立軸を拘束
        end do
! 上記(NSET=FIX2)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX2,1")'
         print'("FIx2,3")'
! 拘束条件の名前を付ける
          print'("*NSET,NSET=FIX3")' !奥の端面の右上の点
! その拘束節点の節点番号
         print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1),','
! 上記(NSET=FIX3)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX3,1,2")'
! 拘束節点の名前を付ける
         print'("*NSET,NSET=FIX4")'
        ! 拘束節点の節点番号
! 奥の端面部分
        do i=2,nfy+1
        ! 奥の端面の縁の1辺を拘束の場合
         print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i,','
        ! 奥の端面の中立軸を拘束の場合
        !print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i-nfx/2,','
        end do
! 上記(NSET=FIX4)の拘束節点の拘束変位
           print'("*BOUNDARY")'
         print'("FIX4,1")'
        print'("*MATERIAL,NAME=EL")'
         print'("*ELASTIC")'
! ヤング率とポアソン比
         !e=210000. ! ヤング率(MPa)
          e=203000.
          e=e*1.d3
         poi=0.3   !  ポアソン比
         print*, e,',',poi
         print'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")'
!
! 載荷条件
! 節点ごとに異なる載荷条件を与える場合をふまえて、
! 載荷条件の種類に対して適当な名前(LASTとか)を付ける
!
!
         print'("*NSET,NSET=LAST")'
      ! 直方体の端面全面に載荷の場合
      !   do i=1,(nfx+1)*(nfy+1)
      !   print'(i6,a,)',(nfx+1)*(nfy+1)*nz+i,','
      !   end do
      !
      !直方体端面の1辺(拘束辺)に載荷の場合
      !   do i=1,nfy+1
      !   print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i,','
      !   end do
      !   
      !直方体端面の中立軸線に載荷の場合
         do i=1,nfy+1
         print'(i6,a,)',(nfx+1)*(nfy+1)*nz+(nfx+1)*i-nfx/2,','
         end do
      !   
      !   
         print'("*STEP")'
         if(nzakutu==0) then
          print'("*STATIC,SOLVER=SPOOLES")'
          print'("*CLOAD")'
          print*,'LAST,3,',(-1.d0/((nfy+1)*(nfx+1)))
         else if(nzakutu==1) then
         print'("*BUCKLE")'
         print*,'5,0.01'
         print'("*CLOAD")'
         print*,'LAST,3,',(-1.d0/((nfy+1)*(nfx+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
