! まずは、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)
! 桁高(y方向)：h(m)
! 梁の軸長(z方向)：ell(m)
!         b=1.2d-2
         b=2.d-2
!         h=2.d-3
         ell=1.d0
! x,y,z方向の要素分割数：nx,ny,nz
         nx=14
         ny=0
         nz=14
!
! 1要素のx,y,z方向の長さ：xn,yn,zn
         xn=b/(2*nx)
         yn=0.000
         zn=ell/nz
!
!         print*,'要素数：',2*nx*nz
!         print*,'節点数：',(2nx+1)*(2nz+1)
! 節点番号ごとの節点座標
! 節点番号：nset
           nset=0
        print*,xn     
           j=0
!        do j=0,2*ny
        do i=0,2*nx
        do k=0,2*nz
           nset=nset+1
           x(nset)=xn*i
           y(nset)=0.000
           z(nset)=zn*k/2
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
!         end do
!
! 要素番号とその要素の6節点の節点番号
!まずは奇数要素から
        nyou1=-1
        print'("*ELEMENT, TYPE=S6, ELSET=Eall")'       
        do k=0,nz-1
        do i=1,nx
         nyou1=nyou1+2
         n1=(2*nx+1)*2*k+2*i-1
         n2=n1+(2*nx+1)
         n3=n2+1
         n4=n1+(2*nx+1)*2
         n5=n3+(2*nx+1)
         n6=n1+(2*nx+1)*2+2
         print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)', &
         nyou1,',',n6,',',n4,',',n1,',',n5,',',n2,',',n3

         end do
         end do
!その後に偶数要素を          
          nyou2=0
          do k=0,nz-1
          do i=1,nx
             nyou2=nyou2+2
           n1=(2*nx+1)*2*k+2*i-1
           n2=n1+1
           n3=n1+2
           n4=n2+(2*nx+1)
           n5=n3+(2*nx+1)
           n6=n1+(2*nx+1)*2+2
!           print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
            print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)', &
           nyou2,',',n3,',',n6,',',n1,',',n5,',',n4,',',n2
!           nyou2,',',n6,',',n3,',',n1,',',n5,',',n2,',',n4
          end do
          end do
          
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の種類に対して適当な名前(FIXとか)を付ける
         print'("*NSET,NSET=FIX")'
! 拘束節点の節点番号
         do k=1,2*nz+1   
         i=(2*nx+1)*(k-1)+1
         print'(i6,a)',i,','
         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")'
         do k=1,2*nz+1
         i=(2*nx+1)*k
         print'(i6,a)',i,','
         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
