! まずは、beam8p.inpと同じデータを、要素番号と節点番号の
! 並び順を変えて任意の要素数で生成するプログラムを作ってみる。
                 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")'
! 桁幅(x方向)：b(m)
! 梁の軸長(z方向)：ell(m)
            b=2.d-1
            ell=1.d0
! x,y,z方向の要素分割数：nx,ny,nz
            nx=14
            ny=0
            nz=8
! 1要素のx,y,z方向の長さ：xn,yn,zn
            xn=b/(2*nx)
            yn=0.000
            zn=ell/(2*nz)
! 節点番号：nset
            nset=0
        do k=0,nx-1
        do i=3*k*nz+2*k+1,(3*k+2)*nz+2*k+1
           nset=nset+1
           x(nset)=xn*2*k
           y(nset)=0.000
           z(nset)=zn*(i-(3*k*nz+2*k+1))
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        nset,',',x(nset),',',y(nset),',',z(nset)
        end do
        do j=(3*k+2)*nz+2*(k+1),3*(k+1)*nz+2*(k+1)
           nset=nset+1 
           x(nset)=xn*(2*k+1)
           y(nset)=0.000
           z(nset)=2*zn*(j-((3*k+2)*nz+2*k+2))
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        nset,',',x(nset),',',y(nset),',',z(nset)
        end do
        end do
!
!
        do i=3*nx*nz+2*k+1,(3*nx+2)*nz+2*nx+1
           nset=nset+1
           x(nset)=xn*2*nx
           y(nset)=0.000
           z(nset)=zn*(i-(3*nx*nz+2*nx+1))
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        nset,',',x(nset),',',y(nset),',',z(nset)
        end do
!        
! 要素番号とその要素の8節点の節点番号
      print'("*ELEMENT, TYPE=S8, ELSET=Eall")'
       do k=1,nx
       do i=1,nz
        n1=3*k*nz+2*i+2*k-1
        n2=3*k*nz+2*i+2*k+1
        n3=(3*k-3)*nz+2*i+2*k-1
        n4=(3*k-3)*nz+2*i+2*k-3
        n5=3*k*nz+2*i+2*k
        n6=(3*k-1)*nz+i+2*k
        n7=(3*k-3)*nz+2*i+2*k-2
        n8=(3*k-1)*nz+i+2*k-1
        print*,(k-1)*nz+i,",",n1,",",n2,",",n3,",",n4,",",n5,",",n6,",",n7,",",n8,","
       end do
       end do
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の種類に対して適当な名前(FIXとか)を付ける
       print'("*NSET,NSET=FIX")'
! 拘束節点の節点番号
       do k=0,nx
       i1=3*k*nz+2*k+1
       print *,i1,","
       end do
       do k=0,nx-1
       i2=(3*k+2)*nz+2*(k+1)
       print *,i2,","
       end do
!
!
! 上記(NSET=FIX)の拘束節点の拘束変位
         print'("*BOUNDARY")'
         print'("FIX,1,6")'
! 節点ごとに異なる拘束条件を与える場合は,
! FIX2とかとして上記を繰り返すのかな(? 調べて)。
         print'("*MATERIAL,NAME=EL")'
         print'("*ELASTIC")'
! ヤング率(MPa)とポアソン比(てことは異方性のときはどうすんのかな)
         print'("70.d3, 0.33")'
         print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")'
!厚さはここで入力
         print'("4.d-3")'
!
! 載荷条件
! 節点ごとに異なる載荷条件を与える場合をふまえて、
! 載荷条件の種類に対して適当な名前(LASTとか)を付ける
!
!
         print'("*NSET,NSET=LAST")'
         i4=2*nz+1
         i5=(3*nx+2)*nz+2*nx+1
         print *,i4,","
         print *,i5,","
!
         do k=1,nx
         i6=3*k*nz+2*k
         print *,i6,","
         end do
!
         do k=0,nx-2
         i7=(3*(k+1)+2)*nz+2*(k+1)+1
         print *,i7,","
         end do
!
!
            print'("*STEP")'
            print'("*STATIC,SOLVER=SPOOLES")'
            print'("*CLOAD")'
!自由端(つまり上で定義したLAST)の全節点(nx1ny1個)の
!y(つまり2)方向に大きさ1の荷重を分配して載荷したいなら、
            write(*,*) 'LAST,1,', 1.d-5/(2*nx+1)
!節点ごとに異なる載荷条件を与える場合は、
!LAST2とかとして上記を繰り返すのかな(?調べて)。
!
             print'("*NODE PRINT,NSET=LAST")'
             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



