      program zgousei1 ver1
      implicit real*8(a-h,o-z)
      dimension nofn(8),b(3,3),e66(6,6),area(100,100)
!
!######## 各断面の諸定数入力 ########
!       write(*,*)'x軸(断面の横の)の長さを入れてください'
!       read(*,*)zl
       zl=6.d-2
!
!       write(*,*)'y軸(断面の縦の)の長さを入れてください'
!       read(*,*)zm
       zm=12.d-2
!
!       write(*,*)'z軸(軸長)の長さを入れてください'
!       read(*,*)ziku
       ziku=1.08d0
!
!       write(*,*)'鋼板のx軸の長さを入れてください'
!       read(*,*)zlk
       zlk=9.d-3
!
!       write(*,*)'上か下の鋼板一本分のy軸の長さを入れてください'
!        write(*,*)'鋼板は上下に入っていることに注意して下さい'
!       read(*,*)zmk
       zmk=4.d-2
!######## 要素数入力 ########
!
!       write(*,*)'x軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!       read(*,*) nx
        nx=8
!
!       write(*,*)'y軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!       read(*,*) ny
        ny=10
!
!       write(*,*)'z軸方向の全要素数を入れて下さい(鋼＋木)'
!       read(*,*) nz
        nz=1
!
!       write(*,*)'鋼板だけのx軸方向の要素数(偶数)を入れてください'
!       read(*,*) nxk
       nxk=2
!       write(*,*)'鋼板だけのy軸方向の要素数(偶数)を入れてください'
!       read(*,*) nyk
       nyk=4
!
!
!######## 出力 #########
        open(7,file='zgousei1.inp')
!
      call zahyou(nx,ny,nz,nxk,nyk,zl,zm,ziku,zmk,zlk,
     & nnofn,zax,zay,zaz)
!
      call youso(nx,ny,nz,nxk,nyk,nofn,nofe)
!
ccccccccccccccccccccccccccc
c 逆行列から9成分を求める c
ccccccccccccccccccccccccccc
c
cccccccccc材料定数ccccccccccccc
ccccccccc 木材の材料定数 cccccccccccccc
c 
c　ヤング率：exx,eyy,ezz
c      write(*,*)'x軸方向のヤング率を入れて下さい'
c      read(*,*) exx
       exx=4.d2
c      write(*,*)'y軸方向のヤング率を入れて下さい'
c      read(*,*) eyy
       eyy=4.d2
c      write(*,*)'z軸方向のヤング率を入れて下さい'
c      read(*,*) ezz
       ezz=1.d4
c
c　ポアソン比：poixy,poixz,poiyz
       poixy=0.016d0
       poixz=poixy
       poiyx=poixy
       poiyz=poixy
       poizx=0.4d0
       poizy=poizx
c
c　せん断弾性係数：gxy,gxz,gxy
       gxy=7.d2
       gxz=gxy
       gyz=gxy
c
c応力−ひずみ行列：e66
          do 10 i=1,6
          do 10 j=1,6
             e66(i,j)=0.d0
10        continue
c  
       e66(1,1)=1.d0/exx
       e66(1,2)=-poixy/exx
       e66(1,3)=-poixz/exx
       e66(2,1)=-poiyx/eyy
       e66(2,2)=1.d0/eyy
       e66(2,3)=-poiyz/eyy
       e66(3,1)=-poizx/ezz
       e66(3,2)=-poizy/ezz
       e66(3,3)=1.d0/ezz
c
c
c
c
c 逆行列を取って、応力−ひずみ行列に変換する
       det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)
     &     +e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)
     &     -e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
       b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
       b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
       b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
       b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
       b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
       b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
       b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
       b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
       b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
             do i=1,3
             do j=1,3
              e66(i,j)=b(j,i)/det
             end do
             end do
c
c
c      write(*,*)'xy面のせん断弾性係数を入れて下さい'
c      read(*,*) gxy
      e66(4,4)=gyz
c      write(*,*)'xz面のせん断弾性係数を入れて下さい'
c      read(*,*) gxz
      e66(5,5)=gxz
c      write(*,*)'yz面のせん断弾性係数を入れて下さい'
c      read(*,*) gyz
      e66(6,6)=gxy
c
c
! ******* 確認のため インプットデッキには書き込まれません ******
!             write(*,*) '応力-ひずみ行列：'
!            do i=1,6
!             write(*,211) (e66(i,j),j=1,6)
!            end do
!211         format(1p6d11.3)
!
cccccccccccccc
c固定端の座標c
cccccccccccccc
       write(7,210)
       do i=1,(nx+1)*(ny+1)
       write(7,*)i,','
       end do
210    format('*NSET,NSET=FIX')
c
cccccccccccccccccccccccccccccccc
cその他境界条件に必要な要素番号c
cccccccccccccccccccccccccccccccc
cccccccc固定端の中立軸の座標cccccccccc
       write(7,520)
       do i =1,nx+1
       write(7,*)ny/2*(nx+1)+i,','
       end do
520    format('*NSET,NSET=CHURITU')
c
cccccccc図心のあるy軸の節点番号cccccccccc
      write(7,'("*NSET,NSET=ZUSIN")')
      do i =1,ny+1 
      write(7,*)(i-1)*(nx+1)+(nx/2)+1,','
      end do
cccccccccc
c境界条件c
cccccccccc
       write(7,220)
       write(7,510)
220    format('*BOUNDARY')
510    format('FIX,3')
c
cccccc固定端の中立軸のyを拘束ccccccc
       write(7,521)
521    format('CHURITU,2')
c
cccccc半分に切断した面のxを拘束ccccccc
c       write(7,531)
c531    format('SETUDAN,1')
c
cccccc図心軸のxを拘束ccccccccccccc
       write(7,'("ZUSIN,1")')
c
!#############
! 材料の設定 #
!#############
!
!#########　鋼材　##########
       write(7,'("*MATERIAL,NAME=STE")') 
!デフォールトは等方性異方性ならTYPE=ANISOをつける
       write(7,'("*ELASTIC")')
!ヤング率、ポワソン比の順
       write(7,'("206000.,.3333")')
!######### 塑性を入れる場合(鋼材) #######
       write(7,'("*PLASTIC")')
       write(7,'("238.,0.")')
       write(7,'("238.,0.5")')
!
!######　直行異方性の場合(集成材)  ######
!
       write(7,231)
231    format('*MATERIAL,NAME=GLU')
       write(7,244)
244    format('*ELASTIC,TYPE=ORTHOTROPIC')
c240    format('*ELASTIC')
       write(7,243) e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','
     & ,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
       write(7,242) e66(4,4)
243    format(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a
     & ,f12.4,a)
242    format(f12.4,',','293.')
c########## 塑性を入れる場合(集成材) #########
!Experimental and analytical study on steel deck-gluam beam
!hybrid bridge behavior : Dr.KISS Lajos
!より集成材の降伏応力を引用
       write(7,'("*PLASTIC")')
       write(7,'("21.696,0.")')
       write(7,'("21.696,0.5")')
c
c
c
!#######################
! 要素に材料特性の付加 #
!#######################
! ここでどの要素がどんな材料か決定
      write(7,'("*SOLID SECTION,ELSET=glulam,MATERIAL=GLU")')
      write(7,'("*SOLID SECTION,ELSET=ukou,MATERIAL=STE")')
      write(7,'("*SOLID SECTION,ELSET=lkou,MATERIAL=STE")')
!
! 
       write(7,'("*STEP")')
       write(7,'("*STATIC,SOLVER=SPOOLES")')
       write(7,'("*CLOAD")')
      call menseki(nx,ny,nz,nxk,zl,zlk,zm,zmk,nyk,
     & nofnp,pakajyu,zakajyu,area)
!
       write(7,'("*NODE PRINT,NSET=Nall")')
       write(7,'("U")')
       write(7,'("*EL PRINT,ELSET=glulam")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=ukou")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=lkou")')
       write(7,'("S")')
       write(7,'("*NODE FILE")')
       write(7,'("U")')
       write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")')
       write(7,'("S")')
       write(7,'("*END STEP")')
!
      close(7)
      end
!###################### ここまでがメイン ###################
!
!######## 座標を出力するサブルーチン######
      subroutine zahyou(nx,ny,nz,nxk,nyk,zl,zm,ziku,zmk,zlk,
     & nnofn,zax,zay,zaz)
      implicit real*8(a-h,o-z)
!
      zxkiten=-1.d0*zl/2.d0
      zykiten=1.d0*zm/2.d0
! xの座標出力
      write(7,'("*NODE,NSET=Nall")')
      nnofn=0
      do k=1,nz+1
      do j=1,ny+1
      do i=1,nx+1
       nnofn=nnofn+1
       nxt=(nx-nxk)/2
       zlt=(zl-zlk)/2.d0
        if(i.le.nxt+1)then
         zax=zxkiten+(i-1)*zlt/nxt
        else if(i.le.nxt+nxk+1)then
         zax=zxkiten+zlt+(i-nxt-1)*zlk/nxk
        else
         zax=zxkiten+zlt+zlk+(i-nx+nxt-1)*zlt/nxt
        end if
!
! yの座標出力
        if(j.le.nyk+1)then
         zay=zykiten+(j-1)*zmk/nyk
        else if(j.le.ny-nyk+1)then
         zay=zykiten+zmk+(j-nyk-1)*(zm-2.d0*zmk)/(ny-2*nyk)
        else
         zay=zykiten+zm-zmk+(j-ny+nyk-1)*zmk/nyk
        end if
!
! zの座標出力
      zaz=(k-1)*ziku/nz
!
      write(7,'(i5,a,F10.7,a,F10.7,a,F10.7)')nnofn,',',zax,',',zay,
     & ',',zaz
      end do
      end do
      end do
!
      return
      end
!
!######## 要素を出力するサブルーチン ########
      subroutine youso(nx,ny,nz,nxk,nyk,nofn,nofe)
      implicit real*8(a-h,o-z)
      dimension nofn(8)
!cccccccccccc 鋼材と木材の節点番号出力 cccccccccccccccc
      do k=1,nz
      do j=1,ny
      do i=1,nx
!****** 上部鋼出力 ******
      nxt=(nx-nxk)/2
      if((i.gt.nxt).and.(i.le.nxt+nxk).and.(j.le.nyk))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=ukou")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
!****** 下部鋼出力 ******
       else if((i.gt.nxt).and.(i.le.nxt+nxk).and.(j.ge.ny-nyk+1))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=lkou")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!****** 集成材出力 ******
       else
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=glulam")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
       end if
      end do
      end do
      end do
!
      return
      end
!
!######## 各要素の面積に応じて載荷するサブルーチン ########
      subroutine menseki(nx,ny,nz,nxk,zl,zlk,zm,zmk,nyk,
     &nofnp,pakjyu,zakajyu,area)
      implicit real*8(a-h,o-z)
      dimension area(100,100),yoko(100),tate(100),ax(100),ay(100)
!******* 各要素の面積を計算 ******
      do j=2,ny+1
      do i=2,nx+1
! xの座標出力
       nxt=(nx-nxk)/2
       zlt=(zl-zlk)/2.d0
        if(i.le.nxt+1)then
         yoko(i)=(i-1)*zlt/nxt
        else if(i.le.nxt+nxk+1)then
         yoko(i)=zlt+(i-nxt-1)*zlk/nxk
        else
         yoko(i)=zlt+zlk+(i-nx+nxt-1)*zlt/nxt
        end if
! yの座標出力
        if(j.le.nyk+1)then
         tate(j)=(j-1)*zmk/nyk
        else if(j.le.ny-nyk+1)then
         tate(j)=zmk+(j-nyk-1)*(zm-2.d0*zmk)/(ny-2*nyk)
        else
         tate(j)=zm-zmk+(j-ny+nyk-1)*zmk/nyk
        end if
      end do
      end do
! 各要素の横の長さ
      yoko(1)=0
      do i=1,nx
      ax(i)=yoko(i+1)-yoko(i)
      end do
! 各要素の縦の長さ
      tate(1)=0
      do j=1,ny
      ay(j)=tate(j+1)-tate(j)
      end do
! 各要素の面積
      do j=1,ny
      do i=1,nx
      area(j,i)=ay(j)*ax(i)
      end do
      end do
!****** 自由端の各節点に載荷 ******
      write(*,'("載荷する荷重を入力して下さい")')
      read*,zkajyu
      zakajyu=zkajyu/zl/zm
!++++++ 4要素と接する節点 ++++++
! 
      do j=1,ny
      do i=1,nx
       pakajyu=(area(j-1,i-1)+area(j-1,i)+area(j,i-1)+area(j,i))/4.d0
     &*zakajyu
       if((i.gt.1).and.(i.lt.nx+1).and.(j.gt.1).and.(j.lt.ny+1))then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
      end if
      end do
      end do
!
!++++++ 2要素と接する節点 ++++++
! 上
      do i =1,nx
       pakajyu=(area(1,i-1)+area(1,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 下
      do i =1,nx
       pakajyu=(area(ny,i-1)+area(ny,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+ny*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 右
      do j =1,ny
       pakajyu=(area(j-1,nx)+area(j,nx))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+1
      write(7,*)nofnp,',','2',',',pakajyu
      end if
      end do
! 左
      do j =1,ny
       pakajyu=(area(j-1,1)+area(j,1))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+nx+1
      write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
!++++++ どことも接しない節点 ++++++
!
      write(7,*)(nx+1)*(ny+1)*nz+1,',','2',',',area(1,1)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+nx+1,',','2',',',
     &area(1,nx)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+ny*(nx+1)+1,',','2',',',
     &area(ny,1)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+ny*(nx+1)+nx+1,',','2',',',
     &area(ny,nx)/4.d0*zakajyu
!
      return
      end
! 
