       program ckhkata ver6
       implicit real*8(a-h,o-z)
       dimension nofn(8),b(3,3),e66(6,6)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  直行異方性材料を解くプログラムだけど等方性の鋼材を解く c
c  ようにしてある。                                       c
c  全解析のデータ入力で勝手に半解析してくれるプログラム。 c
c  荷重条件とか要素とかは全解析のままの数字をどうぞ。荷   c
c  重条件はp,2p,4pを使うタイプ。片持ち梁用。              c
c  この文章をくくっているcのずれは暇なときに直します。    c
c  ver.2006/06/05=ver.6                                   c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
      print'("x軸(断面の横の)の長さを入れてください")'
      read*,zl
      zl=zl/2.d0
!
      print'("y軸(断面の縦の)の長さを入れてください")'
      read*,zm
!       zm=2.d-2
!
      print'("z軸(軸長)の長さを入れてください")'
      read*,ziku
!      ziku=2.d-1
!      
!
!
      print'("x軸方向の要素数(4の倍数以外の偶数)を入れて下さい")'
      read*,nx
      nx=nx/2
!      
      print'("y軸方向の要素数(偶数)を入れて下さい")'
      read*,ny
!
      print'("z軸方向の要素数を入れて下さい")'
      read*,nz
!
!
ccccccccccccccccccccccccccc
c 逆行列から9成分を求める c
ccccccccccccccccccccccccccc
c
cccccccccc材料定数ccccccccccccc
c
c　ヤング率：exx,eyy,ezz
c      write(*,*)'x軸方向のヤング率を入れて下さい'
c      read(*,*) exx
       exx=2.d5
c       exx=1.d4
c      write(*,*)'y軸方向のヤング率を入れて下さい'
c      read(*,*) eyy
       eyy=2.d5
       eyy=exx
c      write(*,*)'z軸方向のヤング率を入れて下さい'
c      read(*,*) ezz
       ezz=2.d5
       ezz=exx
c
c　ポアソン比：poixy,poixz,poiyz
       poixy=2.d-1
       poixy=1.d0/3.d0
       poixz=2.5d-2
       poixz=1.d0/3.d0
       poiyz=3.d-2
       poiyz=1.d0/3.d0
c
c　せん断弾性係数：gxy,gxz,gxy
       gxy=7.5d4
       gxz=gxy
       gyz=gxy
c
c応力−ひずみ行列：e66
          do i=1,6
          do j=1,6
             e66(i,j)=0.d0
          end do
          end do
c  
       e66(1,1)=1.d0/exx
       e66(1,2)=-poixy/exx
       e66(1,3)=-poixz/exx
       e66(2,1)=e66(1,2)
       e66(2,2)=1.d0/eyy
       e66(2,3)=-poiyz/eyy
       e66(3,1)=e66(1,3)
       e66(3,2)=e66(2,3)
       e66(3,3)=1.d0/ezz
c
c
c
c ヤング率やポアソン比を書き写してから
      exx=1.d0/e66(1,1)
      eyy=1.d0/e66(2,2)
      ezz=1.d0/e66(3,3)
      poixy=-e66(1,2)*exx
      poixz=-e66(1,3)*exx
      poiyz=-e66(2,3)*eyy
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)=gxy
c      write(*,*)'xz面のせん断弾性係数を入れて下さい'
c      read(*,*) gxz
      e66(5,5)=gxz
c      write(*,*)'yz面のせん断弾性係数を入れて下さい'
c      read(*,*) gyz
      e66(6,6)=gyz
c
c
       open(7,file='ckhkata.inp')
cccccccccc
c座標出力c
cccccccccc
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 確認のため
             write(*,*) '応力-ひずみ行列：'
            do i=1,6
             write(*,211) (e66(i,j),j=1,6)
            end do
211         format(1p6d11.3)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
       write(7,200)
       nnofe=0
       do k=1,nz+1
       do j=1,ny+1
       do i=1,nx+1
        nnofe=nnofe+1
        zax=(i-1)*zl/nx
        zay=(j-1)*zm/ny
        zaz=(k-1)*ziku/nz
       write(7,100)nnofe,',',zax,',',zay,',',zaz
       end do
       end do
       end do
200    format('*NODE,NSET=Nall')
100    format(i5,a,F6.4,a,F6.4,a,F6.4)
c  別に書かなくてもいいのだろうが、format文の書き方は
c  http://www.g.dendai.ac.jp/lecture/fortran/fort04.htmlを見て勉強。
c  探した中で は一番分かりやすい。
cccccccccccccc
c節点番号出力c
cccccccccccccc
! nofnは number of node の略。ノフンて嫌な響だけど...
       do k=1,nz
       do j=1,ny
       do i=1,nx
        node=(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,110)
        write(7,*)node,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
       end do
       end do
       end do
110    format('*ELEMENT,TYPE=C3D8,ELSET=Eall')
c
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,','
c       write(7,*)i,','
       end do
520    format('*NSET,NSET=CHURITU')
c
cccccccc切断面の要素番号cccccccccccccc
       write(7,530)
       do k=1,nz+1
       do j=1,ny+1
       write(7,*) (nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+1,','
       end do
       end do
530    format('*NSET,NSET=SETUDAN')
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
       write(7,531)
531    format('SETUDAN,1')
c
cccccc
c材料c
cccccc
       write(7,230)
230    format('*MATERIAL,NAME=EL')
cccccccccc
c弾性特性c
cccccccccc
       write(7,'("*ELASTIC")')
cデフォールトは等方性.異方性ならTYPE=ANISOをつける
       write(7,'(("207000.,"),f6.5)')1.d0/3.d0
!       write(7,250)
c250    format('210000.,.3')
!250    format('200000.,.3333')
c
c
c
c
ccccccc　直行異方性の場合  ccccccc
c       write(7,240)
c240    format('*ELASTIC,TYPE=ORTHOTROPIC')
c       write(7,243) e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','
c     & ,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5)
c       write(7,242) e66(4,4)
c243    format(f12.1,a,f12.1,a,f12.1,a,f12.1,a,f12.1,a,f12.1,a,f12.1,a
c     & ,f12.1)
c242    format(f12.1,',','294')
c
c
c
c
cccccccccccccccccccccccc
c 材料特性を節点に付加 c
cccccccccccccccccccccccc
       write(7,260)
260    format('*SOLID SECTION,ELSET=Eall,MATERIAL=EL')
cccccccccccccc
c自由端の座標c
cccccccccccccc
       write(7,270)
       nnofe=(nx+1)*(ny+1)*(nz+1)-(nx+1)*(ny+1)
       do i=1,(nx+1)*(ny+1)
       nnofe=nnofe+1
       write(7,*)nnofe,','
       end do
270    format('*NSET,NSET=LAST')
c
c
c
c
       write(7,300)
300    format('*STEP')
       write(7,310)
310    format('*STATIC,SOLVER=SPOOLES')
       write(7,320)
320    format('*CLOAD')
c
cccccccccc
c載荷条件c
cccccccccc
cccccccccccccccccccccccccccccc
c各節点に同じ荷重をかける場合c
cccccccccccccccccccccccccccccc
c       zkajyuu=9.d0/(2*nx+1)/(ny+1)
c       do i =1,ny+1
c       khaji=(nx+1)*(ny+1)*nz+i*(nx+1)
c       write(7,*)khaji,',','2',',',zkajyuu/2.d0
c       end do
c       do j=1,ny+1
c       do i=1,nx
c       knoko=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i
c       write(7,*)knoko,',','2',',',zkajyuu
c       end do
c       end do
cccccccccccccccccccccccccccccccccc
c各節点に荷重を変えて載荷する場合c
cccccccccccccccccccccccccccccccccc
c
ccccccc p,2p,4pに分けて載荷するけど、そのpの求め方(単位はMN)cccccc
      print*,'載荷する荷重を入力して下さい(MN)'
      read'(f10.6)',zkajyu
c       zkajyu=2.d-3
! 拘束箇所に載荷できないので、一つ置きに載荷
      nzkajyu=nx*ny !nzkajyuは載荷する節点の数
       pkajyu=zkajyu/nzkajyu/2.d0
c
ccccccc p載荷する節点番号と荷重の方向と大きさの出力cccccccccccc
c
       do i=1,ny+1,ny
       write(7,*)(nx+1)*(ny+1)*nz+(i-1)*(nx+1)+nx+1,',',2,',',pkajyu
       end do
c
ccccccc 2p載荷する節点番号と荷重の方向と大きさの出力cccccccccccc
       p2kajyu=pkajyu*2.d0
       do j=1,ny+1,ny
       do i =2,nx,2
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',p2kajyu
       end do
       end do
c
       do i =3,ny,2
       write(7,*) (nx+1)*(ny+1)*nz+(i-1)*(nx+1)+nx+1,',','2',',',p2kajyu
       end do
c
ccccccc 4p載荷する節点番号と荷重の方向と大きさの出力cccccccccccccc
       p4kajyu=pkajyu*4.d0
       do j =3,ny,2
       do i =2,nx,2
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',p4kajyu
       end do
       end do
c
       write(7,'("*NODE PRINT,NSET=Nall")')
       write(7,350)
350    format('U')
       write(7,360)
360    format('*EL PRINT,ELSET=Eall')
       write(7,370)
370    format('S')
       write(7,380)
380    format('*NODE FILE')
       write(7,350)
       write(7,390)
390    format('*EL FILE,POSITION=AVERAGED AT NODES')
       write(7,370)
       write(7,400)
400    format('*END STEP')
c
       close(7)
c
c
c
c
c
       end
