! まずは、beam8p.inpと同じデータを、要素番号と節点番号の
! 並び順を変えて任意の要素数で生成するプログラムを作ってみる。
                 implicit real*8(a-h,o-z)
                 dimension x1(100000),x2(100000),x3(100000),x4(100000),y(100000),z(100000)
! 節点座標：x(節点番号),y(節点番号),z(節点番号)
           print'("**")'
           print'("**beam8p.inpと同じデータを要素番号と節点番号の")'
           print'("**並び順(と要素数)を変えたデータ")'
           print'("**")'
           print'("*NODE,NSET=Nall")'
!z方向の長さを入力
            d=2.d-1
! x,y,z方向の要素分割数：nx,ny,nz
!辺の長さ
           b=3.0
!分割数を代入する
!偶数でないとダメ
           nx=4
           ny=20
!            nz=20
!分割数で割った値が6角形一つの長さとなる
! 1要素のx,y,z方向の長さ：xn,yn,zn
       open(7,file='yousosuu.txt')
       write(7,*) 4*nx+4*nx*(ny-1)+2*nx+2*nx+4*nx*ny
       write(7,*) 2*nx,2*ny
       close(7)
!
!dは厚さ、xnはハニカムを構成する六角形の一辺だが
!yousosuu2.txtの中にx方向の長さとy方向の長さを
!自動的に出力するので任意の長さに自分で決める。
            xn=b/nx
       open(7,file='yousosuu2.txt')
       write(7,*) ,"x方向の長さは",2*xn*nx+xn*(nx-1)
       write(7,*) ,"y方向の長さは",xn*3**(0.5)*ny
       close(7)
!            yn=0.000
!            zn=ell/(2*nz)
! 節点番号
!一列目
            nset=0
        do k=1,ny+1
        do j=1,2
        do i=1,nx
           nset=nset+1
!           nset1=nset+2
          x1(nset)=3*xn*i-3*xn
          x2(nset)=3*xn*i-2*xn
           y(nset)=(xn*(3)**(0.5))*(k-1)
           z(nset)=d*(j-1)
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
!         nset,',',x1(nset),',',y(nset),',',z(nset) 
        (2*i-1)+2*nx*(j-1)+4*nx*(k-1),',',x1(nset),',',y(nset),',',z(nset)
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
!         nset,',',x2(nset),',',y(nset),',',z(nset)
        2*i+2*nx*(j-1)+4*nx*(k-1),',',x2(nset),',',y(nset),',',z(nset)
        end do
        end do
        end do
!2列目
            nset=0
        do k=1,ny
        do j=1,2
        do i=1,nx
           nset=nset+1
          x3(nset)=3*xn*i-3*xn-xn*(0.5)
          x4(nset)=3*xn*i-2*xn+xn*(0.5)
!          x4(nset)=3*xn*i-(2*xn-xn*(0.5))
!           x(nset)=(r*3/2)-((r*1/2)*(-1)**i)
           y(nset)=xn*(3)**(0.5)/2+(xn*(3)**(0.5))*(k-1)
           z(nset)=d*(j-1)
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        4*nx-1+4*nx*(ny-1)+(2*i-1)+2*nx*(j-1)+4*nx*k+1,',',x3(nset),',',y(nset),',',z(nset)
        print'(i6,a,f9.6,a,f9.6,a,f9.6)', &
        4*nx+4*nx*(ny-1)+2*i+2*nx*(j-1)+4*nx*k,',',x4(nset),',',y(nset),',',z(nset)
        end do
        end do
        end do
!
!
!---------------------------------------------------------------------------------------------------
!
!
! 要素番号とその要素の8節点の節点番号
        print'("*ELEMENT, TYPE=S3, ELSET=Eall")'
!1
        do k=1,ny/2
        do i=1,nx
        n1=2*i-1+(8*nx)*(k-1)
        n2=n1+1
        n3=n1+nx*2
        n4=n3+1
       print*,2*i-1+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
       print*,2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
        end do
!2
        do i=1,nx
        n1=2*nx*2*(ny+1)+2*i-1+(8*nx)*(k-1)
        n2=2*i-1+(8*nx)*(k-1)
        n3=n1+nx*2
        n4=n2+2*nx
       print*,nx*2+2*i-1+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
       print*,nx*2+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
       end do
!3
       do i=1,nx
       n1=2*i+(8*nx)*(k-1)
       n2=2*nx*2*(ny+1)+2*i-1+1+(8*nx)*(k-1)
       n3=2*nx+2*i-1+1+(8*nx)*(k-1)
       n4=n2+nx*2
      print*,nx*2+2*nx-1+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
      print*,nx*2+2*nx+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
       end do
!4
      do i=1,nx
       n1=2*nx*2+2*i-1+(8*nx)*(k-1)
       n2=2*nx*2*(ny+1)+2*i-1+(8*nx)*(k-1)
       n3=n1+2*nx
       n4=n2+2*nx
      print*,nx*2+4*nx-1+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
      print*,nx*2+4*nx+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
      end do
!5
     do i=1,nx
     n1=2*nx*2+2*i-1+(8*nx)*(k-1)
     n2=n1+1
     n3=n1+nx*2
     n4=n3+1
     print*,nx*2+6*nx-1+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
     print*,nx*2+6*nx+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
     end do
!6
    do i=1,nx
    n1=2*nx*2+2*i-1+1+(8*nx)*(k-1)
    n2=2*nx*2*(ny+1)+2*i-1+1+(8*nx)*(k-1)
    n3=n1+nx*2
    n4=n2+nx*2
    print*,nx*2+8*nx-1+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
    print*,nx*2+8*nx+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
    end do
!7
    do i=1,nx-1
    n1=2*nx*2*(ny+1)+2*i-1+1+(8*nx)*(k-1)
    n2=n1+1
    n3=n1+nx*2
    n4=n3+1
    print*,nx*2+10*nx-1+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
    print*,nx*2+10*nx+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
    end do
!8
    do i=1,nx
    n1=2*nx*2*(ny+1)+2*i-1+1+nx*4-1+(8*nx)*(k-1)
    n2=2*nx*2+2*i-1+(8*nx)*(k-1)
    n3=n1+nx*2
    n4=n2+nx*2
    print*,nx*2+12*nx-3+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
    print*,nx*2+12*nx-2+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
    end do
!9
   do i=1,nx
   n1=2*nx*2+2*i-1+1+(8*nx)*(k-1)
   n2=2*nx*2*(ny+1)+2*i-1+1+nx*4-1+1+(8*nx)*(k-1)
   n3=n1+nx*2
   n4=n2+nx*2
   print*,nx*2+14*nx-3+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
   print*,nx*2+14*nx-2+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
   end do
!10
   do i=1,nx-1
   n1=2*nx*2*(ny+1)+2*i-1+1+nx*4-1+1+(8*nx)*(k-1)
   n2=n1+1
   n3=n1+nx*2
   n4=n3+1
   print*,nx*2+16*nx-3+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
   print*,nx*2+16*nx-2+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
   end do
!11
   do i=1,nx
   n1=(4*nx+2*i-1)*2-(2*i-1)+(8*nx)*(k-1)
   n2=2*nx*2*(ny+1)+2*i-1+1+nx*4-1+(8*nx)*(k-1)
   n3=n1+nx*2
   n4=n2+nx*2
   print*,nx*2+18*nx-5+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
   print*,nx*2+18*nx-4+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
   end do
!12
   do i=1,nx
   n1=(4*nx+2*i-1)*2-(2*i-1)+1+(8*nx)*(k-1)
   n2=2*nx*2*(ny+1)+2*i-1+1+nx*4-1+1+(8*nx)*(k-1)
   n3=n1+nx*2
   n4=n2+nx*2
   print*,nx*2+20*nx-5+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n1,",",n2,",",n3
   print*,nx*2+20*nx-4+2*i+(nx*2+20*nx-5+2*nx+1)*(k-1),",",n4,",",n3,",",n2
   end do
   end do
!最後はこれで完成
!13
        do i=1,nx
        n1=2*i-1+(8*nx)*(ny/2)
        n2=n1+1
        n3=n1+nx*2
        n4=n3+1
       print*,nx*2+20*nx-5+2*nx+(nx*2+20*nx-5+2*nx+1)*(ny/2-1)+2*i,",",n1,",",n2,",",n3
       print*,nx*2+20*nx-4+2*nx+(nx*2+20*nx-5+2*nx+1)*(ny/2-1)+2*i,",",n4,",",n3,",",n2
        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
