      program cghkata ver a2!(2006/11/28)
      implicit real*8(a-h,o-z)
      dimension nofn(8),b(3,3),e66(6,6),area(100,100),e77(6,6)
!
!######## 各断面の諸定数入力 ########
!       print*,'x軸(断面の横の)の長さを入れてください'
!       read*,zl
       zl=9.6*2d-2
       zl=zl/2.d0
!
!       print*,'y軸(断面の縦の)の長さを入れてください'
!       read*,zm
       zm=9.6d-2
!
!       print*,'z軸(軸長)の長さを入れてください'
!       read*,ziku
       ziku=ziku/2.d0
       ziku=1.0*2d0
!
!######## 要素数入力 ########
!
!       print*,'x軸方向の全要素数(4以上の偶数)を入れて下さい'
!       read*,nx
       nx=4
       nx=nx/2
!
!       print*,'y軸方向の全要素数(偶数)を入れて下さい'
!       read*, ny
       ny=6
!
!       print*,'z軸方向の全要素数を入れて下さい'
!       read*,nz
       nz=460
!
!
!######## 出力 #########
        open(7,file='cghkata.inp')
!
!座標読み込み
      call zahyou(nx,ny,nz,zl,zm,ziku,nnofn,zax,zay,zaz)
!
!要素番号読み込み
      call youso(nx,ny,nz,nofn,nofe)
!
!######## 境界条件 ########
!境界条件に必要な座標
!
!固定端の座標
       write(7,'("*NSET,NSET=FIX")')
       do i=1,(nx+1)*(ny+1)
       write(7,*)i,','
       end do
!固定端の中立軸の座標
       write(7,'("*NSET,NSET=CHURITU")')
       do i =1,nx+1
       write(7,*)ny/2*(nx+1)+i,','
       end do
!図心を含むようにyz平面で切断した面の座標
      write(7,'("*NSET,NSET=SETUDAN")')
      do j=1,nz+1
      do i =1,ny+1 
      write(7,*)(j-1)*(ny+1)*(nx+1)+(i-1)*(nx+1)+1,','
      end do
      end do
!
      write(7,'("*BOUNDARY")')
      write(7,'("FIX,3")')
      write(7,'("CHURITU,2")')
      write(7,'("SETUDAN,1")')
!
!######## 材料の設定 ########
!******** 集成材の材料設定(圧縮) ********
      write(7,'("*MATERIAL,NAME=GLUP")')
      write(7,'("*ELASTIC,TYPE=ORTHOTROPIC")')
! 
!　ヤング率：exx,eyy,ezz
!      print*,'z軸方向のヤング率を入れて下さい'
!      read*, ezz
       ezz=9.00d3
!      print*,'x軸方向のヤング率を入れて下さい'
!      read*, exx
       exx=ezz/25.d0
!      print*,'y軸方向のヤング率を入れて下さい'
!      read*, eyy
       eyy=ezz/25.d0
!
!　ポアソン比：poixy,poixz,poiyz
       poixy=0.016d0
       poixz=poixy
       poiyx=poixy
       poiyz=poixy
       poizx=0.4d0
       poizy=poizx
!
!　せん断弾性係数：gxy,gxz,gxy
       gxy=6.00d2
       gxz=gxy
       gyz=gxy
!
      call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,
     & gxz,gxy,poiyx,poizx,poizy)
       write(7,'(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,
     & a ,f12.4,a)')
     &  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,'(f12.4,a,a)')e66(4,4),',','294'
!++++++++++ 塑性を入れる場合(集成材・圧縮) +++++++++
!Experimental and analytical study on steel deck-gluam beam
!より集成材の降伏応力を引用
       write(7,'("*PLASTIC")')
       write(7,'("21.696,0.")')
       write(7,'("21.696,0.5")')
!******** 集成材の材料設定(引張) ********
       write(7,'("*MATERIAL,NAME=GLUT")')
       write(7,'("*ELASTIC,TYPE=ORTHOTROPIC")')
       write(7,'(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a
     & ,f12.4,a)')
     &  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,'(f12.4,a,a)')e66(4,4),',','294'
!
!#######################
! 要素に材料特性の付加 #
!#######################
! ここでどの要素がどんな材料か決定
      write(7,'("*SOLID SECTION,ELSET=glulamp,MATERIAL=GLUP")')
      write(7,'("*SOLID SECTION,ELSET=glulamt,MATERIAL=GLUT")')
!
! 
       write(7,'("*STEP")')
       write(7,'("*STATIC,SOLVER=SPOOLES")')
       write(7,'("*CLOAD")')
!載荷条件読み込み
      call menseki(nx,ny,nz,zl,zm,
     & nofnp,pakajyu,zakajyu,area)
!
       write(7,'("*NODE PRINT,NSET=Nall")')
       write(7,'("U")')
       write(7,'("*EL PRINT,ELSET=glulamp")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=glulamt")')
       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,zl,zm,ziku,nnofn,zax,zay,zaz)
      implicit real*8(a-h,o-z)
!
      zykiten=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
        if(i.le.2)then
         zax=(i-1)*zl/(2*nx-1)
        else
         zax=zl/(2*nx-1)+(i-2)*zl*2.d0/(2*nx-1)
        end if
!
! yの座標出力
      zay=zykiten+(j-1)*zm/ny
!
! 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,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
        if(j.le.ny/2)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=glulamp")')
      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=glulamt")')
      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,zl,zm,
     &nofnp,pakajyu,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の座標出力
        if(i.le.2)then
         yoko(i)=(i-1)*zl/(2*nx-1)
        else
         yoko(i)=zl/(2*nx-1)+(i-2)*zl*2.d0/(2*nx-1)
        end if
!
! yの座標出力
         tate(j)=(j-1)*zm/ny
      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
      zkajyu=1.0d-3
      zakajyu=zkajyu/zl/zm
!片持ち梁だから2で割る、さらに半解析だから2で割る
      zakajyu=zakajyu/2.d0/2.d0*-1.d0
!++++++ 4要素と接する節点 ++++++
! 対称面を含む要素の節点
      do j=2,ny
       pakajyu=(2.d0*area(j-1,1)+area(j-1,2)+2.d0*area(j,1)+area(j,2))
     &/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+2
       write(7,*)nofnp,',','2',',',pakajyu
      end do
! その他 
      do j=2,ny
      do i=3,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要素と接する節点 ++++++
! 上(対称面を含む要素)
       pakajyu=(2.d0*area(1,1)+area(1,2))/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nz+2
       write(7,*)nofnp,',','2',',',pakajyu
! 上
      do i =3,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
! 下(対称面を含む要素)
       pakajyu=(2.d0*area(ny,1)+area(ny,2))/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nz+ny*(nx+1)+2
       write(7,*)nofnp,',','2',',',pakajyu
! 下
      do i =3,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 =2,ny
       pakajyu=(area(j-1,nx)+area(j,nx))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+j*(nx+1)
      write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
!++++++ どことも接しない節点 ++++++
!
      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)+nx+1,',','2',',',
     &area(ny,nx)/4.d0*zakajyu
!
      return
      end
! 
!
!######## 直交異方性の定数を出すプログラム ########
      subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,
     & gxz,gxy,poiyx,poizx,poizy)
      implicit real*8(a-h,o-z)
      dimension e66(6,6),b(3,3)
!応力−ひずみ行列：e66
          do i=1,6
          do j=1,6
             e66(i,j)=0.d0
          end do
          end do
!  
       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
!
!****** 逆行列を取って、応力−ひずみ行列に変換する******
       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
!
!****** せん断補正係数を入力 ******
!      print*,'xy面のせん断弾性係数を入れて下さい'
!      read*, gxy
      e66(4,4)=gyz
!      print*,'xz面のせん断弾性係数を入れて下さい'
!      read*, gxz
      e66(5,5)=gxz
!      print*,'yz面のせん断弾性係数を入れて下さい'
!      read*, gyz
      e66(6,6)=gxy
!
!
      return
      end
