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