       program C3D8 hasira ver4
       implicit real*8(a-h,o-z)
       dimension nofe(8),b(3,3),e66(6,6),c(3,3)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  全解析のデータ入力で勝手に1/8解析してくれるプログラム。        c
c  荷重条件とか要素とかは全解析のままの数字をどうぞ。荷           c
c  重条件はp,2p,4pを使うタイプ。柱用。引張試験用?               c
c  この文章をくくっているcのずれは暇なときに直します。            c
c  ver.2005/12/28=ver.4                               c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c       write(*,*)'x軸(断面の横の)の長さを入れてください'
c       read(*,*)zl
       zl=2.d-2
       zl=zl/2.d0
c       write(*,*)'y軸(断面の縦の)の長さを入れてください'
c       read(*,*)zm
       zm=2.d-2
       zm=zm/2.d0
c       zm=5.d0
c       write(*,*)'z軸(軸長)の長さを入れてください'
c       read(*,*)ziku
       ziku=2.d-1
       ziku=ziku/2.d0
c       
c
c
c       write(*,*)'x軸方向の要素数(偶数)を入れて下さい'
c       read(*,*) nx
       nx=2
       nx=nx/2
c       write(*,*)'y軸方向の要素数(偶数)を入れて下さい'
c       read(*,*) ny
       ny=2
       ny=ny/2
c       write(*,*)'z軸方向の要素数を入れて下さい'
c       read(*,*) nz
       nz=2
       nz=nz/2
c       write(*,*)どの方向を拘束しますか
c
       nzen=nx*ny*nz
c       write(*,*)'ヤング率をPaで入れてください'
c       read(*,*) nyang
c       write(*,*)'ポワソン比を小数点以下の数字だけ入れてください'
c       write(*,*)'例えば0.3なら3'
c       read(*,*) npoi
c
ccccccccccccccccccccccccccc
c 逆行列から9成分を求める c
ccccccccccccccccccccccccccc
c
cccccccccc材料定数ccccccccccccc
c
c　ヤング率：exx,eyy,ezz
c      write(*,*)'x軸方向のヤング率を入れて下さい'
c      read(*,*) exx
       exx=5.d2
       exx=2.d5
c       exx=1.d4
c      write(*,*)'y軸方向のヤング率を入れて下さい'
c      read(*,*) eyy
       eyy=1.d3
       eyy=2.d5
       eyy=exx
c      write(*,*)'z軸方向のヤング率を入れて下さい'
c      read(*,*) ezz
       ezz=1.d4
       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 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)=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 20 i=1,3
             do 20 j=1,3
              e66(i,j)=b(j,i)/det
20           continue
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='hasira.inp')
cccccccccc
c座標出力c
cccccccccc
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 確認のため
c             write(*,*) '応力-ひずみ行列：'
c            do 211 i=1,6
c             write(*,212) (e66(i,j),j=1,6)
c211         continue
c212         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
       do k=1,nz
       do j=1,ny
       do i=1,nx
        node=(k-1)*nx*ny+nx*(j-1)+i
        nofe(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofe(2)=nofe(1)+1
        nofe(3)=nofe(1)+(nx+1)
        nofe(4)=nofe(3)+1
        nofe(5)=nofe(1)+(nx+1)*(ny+1)
        nofe(6)=nofe(2)+(nx+1)*(ny+1)
        nofe(7)=nofe(3)+(nx+1)*(ny+1)
        nofe(8)=nofe(4)+(nx+1)*(ny+1)
        write(7,110)
        write(7,*)node,',',nofe(1),',',nofe(3),',',nofe(7),
     & ',',nofe(5),',',nofe(2),',',nofe(4),',',nofe(8),',',nofe(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
c       write(7,520)
c       do i =1,nx+1
c       write(7,*)ny/2*(nx+1)+i,','
c       write(7,*)i,','
c       end do
c520    format('*NSET,NSET=CHURITU')
c
cccccccc切断面の要素番号cccccccccccccc
c
cccccccc yz平面 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=SETUYZ')
c
cccccccc xz平面 cccccccccccccc
       write(7,540)
       do k=1,nz+1
       do i=1,nx+1
       write(7,*) (k-1)*(nx+1)*(ny+1)+i,','
       end do
       end do
540    format('*NSET,NSET=SETUXZ')
cccccccccc
c境界条件c
cccccccccc
       write(7,220)
       write(7,510)
220    format('*BOUNDARY')
510    format('FIX,3')
c
cccccc固定端の中立軸のyを拘束ccccccc
c       write(7,521)
c521    format('CHURITU,2')
c
cccccc 半分に切断したyz面のxを拘束 ccccccc
       write(7,531)
531    format('SETUYZ,1')
c
cccccc 半分に切断したxz面のyを拘束 ccccccc
       write(7,541)
541    format('SETUXZ,2')
cccccc
c材料c
cccccc
       write(7,230)
230    format('*MATERIAL,NAME=EL')
cccccccccc
c弾性特性c
cccccccccc
       write(7,240)
240    format('*ELASTIC')
cデフォールトは等方性異方性ならTYPE=ANISOをつける
       write(7,250)
c250    format('210000.,.3')
250    format('207000,.3333')
c
c250    format('nyang.,.npoi')
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
cccccccc 塑性を入れるとき ccccccccccc
       write(7,244)
c244    format('*PLASTIC')
244    format('*PLASTIC,HARDENING=KINEMATIC')
       write(7,*)'238.,0.'
       write(7,*)'242.,0.01405'
c       write(7,*)'291.,0.010'
       write(7,*)'319.625,0.02905'
c上の定数はKhaledさんの
c Spatial ultimate load analysis of thin-walled and its
c application to the design of channel members より抜粋
c なお、三行目の290.75,0.01は、論文に書かれていなかったため
c 最初に相当塑性ひずみを0.01とおき、Esを乗じて応力を出している
c       write(7,*)'1000.,0.15'
c
c
c
cccccccccccccccccc
c材料特性の3D化？c
cccccccccccccccccc
       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
c       zkajyu=2.d-3
c       zkajyu=2.d-5
       write(*,*)'どのくらいの力で引っ張りますか?単位はMN'
       read(*,*)zkajyu
       pkajyu=zkajyu/4.d0/nx/ny/4.d0
c
ccccccc p載荷する節点番号と荷重の方向と大きさの出力cccccccccccc
       do j =1,ny+1,ny
       do i =1,nx+1,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','3',','
     &,pkajyu
       end do
       end do
c
ccccccc 2p載荷する節点番号と荷重の方向と大きさの出力cccccccccccc
       p2kajyu=pkajyu*2.d0
       do j =1,ny+1,ny
       do i =2,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','3',','
     &,p2kajyu
       end do
       end do
c
       do j =2,ny
       do i =1,nx+1,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','3',','
     & ,p2kajyu
       end do
       end do
c
ccccccc 4p載荷する節点番号と荷重の方向と大きさの出力cccccccccccccc
       p4kajyu=pkajyu*4.d0
       do j =2,ny
       do i =2,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','3',','
     & ,p4kajyu
       end do
       end do
c
c       write(7,330)
c330    format('LAST,2,0.36')
       write(7,340)
340    format('*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
