! 両端支持に任意の面積で荷重をかけ解析するプログラム     
     
      implicit real*8(a-h,o-z)
      dimension x(10000),y(10000),z(10000),e66(6,6)
! コンパイルして実行するとhari.inpができる。      

      open(7,file='harioa.inp')

! 桁幅(x方向)：b(m)
! 桁高(y方向)：h(m)
! 梁の軸長(z方向)：ell(m)
      b=9.6d-2
      h=9.6d-2
      ell=1.0d0

! x,y,z方向の要素分割数：nx,ny,nz
      nx=2
      ny=6
      nz=400



      write(7,'("*NODE, NSET=Nall")')

! 1要素のx,y,z方向の長さ：xn,yn,zn
      xn=b/nx
      yn=h/ny
      zn=ell/nz

! 節点番号ごとの節点座標
! 節点番号：nset
           nset=0
        do k=0,nz
        do j=0,ny
        do i=0,nx
           nset=nset+1
           x(nset)=xn*i
           y(nset)=yn*j
           z(nset)=zn*k
        write(7,'(i6,a,f9.3,a,f9.3,a,f9.3)')
     &         nset,',',x(nset),',',y(nset),',',z(nset)
         end do
         end do
         end do
! 要素番号とその要素の8節点の節点番号
         nyou=0
        do k=1,nz
        do j=1,ny
        do i=1,nx
         nyou=nyou+1
         n1=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
         n2=n1+1
         n3=n1+(nx+1)
         n4=n3+1
         n5=n1+(nx+1)*(ny+1)
         n6=n2+(nx+1)*(ny+1)
         n7=n3+(nx+1)*(ny+1)
         n8=n4+(nx+1)*(ny+1)
         write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Eall")')
         write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')
     &   nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
         end do
         end do
         end do

! 始端の固定支承の節点
         write(7,'("*NSET,NSET=KOTEI")')
! 節点の抜き出し
         nx2ny2=(nx+1)*ny+1
         nx1ny1=(nx+1)*(ny+1)
         do i=nx2ny2,nx1ny1
         write(7,'(i6,a)'),i,','
         end do
  

! 支承の面積を決定する。1大きくなるごとに１列増える。
       nn=5

! 始端のローラー支承の節点 
         write(7,'("*NSET,NSET=SITAN")')
! 節点の抜き出し
         aa=nx2ny2
         bb=nx1ny1
         do l=1,nn       
         aa=aa+nx1ny1
         bb=bb+nx1ny1
         do i=aa,bb
         write(7,'(i6,a)'),i,','
         end do
         end do

! 終端のローラー支承の節点         
        write(7,'("*NSET,NSET=SHUTAN")')
! 節点の抜き出し
         cc=(nx+1)*(ny+1)*(nz-nn)-nx
         dd=(nx+1)*(ny+1)*(nz-nn)
         do l=1,nn       
         cc=cc+nx1ny1
         dd=dd+nx1ny1
         do i=cc,dd
         write(7,'(i6,a)'),i,','
         end do
         end do

! 終端の一番端のローラー支承の節点         
        write(7,'("*NSET,NSET=ROLL")')
! 節点の抜き出し
         ee=(nx+1)*(ny+1)*(nz+1)-nx
         ff=(nx+1)*(ny+1)*(nz+1)
         do i=ee,ff
         write(7,'(i6,a)'),i,','
         end do


! 拘束について
! 1でx軸を、2でy軸を、3でz軸を固定する。
! 4~6でそれぞれの回転も固定できる。
         write(7,'("*BOUNDARY")')
         write(7,'("KOTEI,1,3")')
         write(7,'("SITAN,1,2")')
         write(7,'("SHUTAN,1,2")')
         write(7,'("ROLL,1,2")')

! 材料の設定 

!ここで直交異方性を指定している。
      write(7,'("*MATERIAL,NAME=EL")')
      write(7,'("*ELASTIC,TYPE=ORTHOTROPIC")')

! ヤング率：exx,eyy,ezz

!      z軸方向のヤング率
       ezz=9.00d3

!      x軸方向のヤング率を入れて下さい
       exx=ezz/25.d0
!      y軸方向のヤング率を入れて下さい
       eyy=ezz/25.d0
!
! ポアソン比：poixy,poixz,poiyz
! poixyとpoizxは逆数の関係
       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'
       write(7,'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")')

! 載荷条件
! 載荷する荷重の節点   

       
! 荷重をかける面積を決定する。梁の中央を基準に左右に１列づつ増えていく。
         mm=0

         write(7,'("*NSET,NSET=KAJU")')
         
         aaa=(1+((nx+1)*(ny+1)*nz+1))/2
         bbb=(1+((nx+1)*(ny+1)*nz+1))/2+nx
         ccc=aaa
         ddd=bbb
         eee=aaa
         fff=bbb        
         
         do i=aaa,bbb
         write(7,'(i6,a)'),i,','
         end do
         
         do l=1,mm
         ccc=ccc-nx1ny1
         ddd=ddd-nx1ny1
         do i=ccc,ddd
         write(7,'(i6,a)'),i,','
         end do
         end do
         
         do l=1,mm
         eee=eee+nx1ny1
         fff=fff+nx1ny1        
         do i=eee,fff
         write(7,'(i6,a)'),i,','
         end do
         end do


         write(7,'("*STEP")')
         write(7,'("*STATIC,SOLVER=SPOOLES")')
! CLOADで集中荷重を設定
         write(7,'("*CLOAD")')


!載荷する面積
         menseki=(nx+1)*(2*mm+1)
! 荷重の大きさ
! 載荷する荷重の方向　y軸なので2
! 載荷する荷重の大きさを面積で割ることで、全体に等分布に載荷される。
         write(7,*) 'KAJU,2,', 1.0d-3/menseki



         write(7,'("*NODE PRINT,NSET=Nall")')
         write(7,'("U")')
         write(7,'("*EL PRINT,ELSET=Eall")')
         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")')


         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
!
!****** せん断補正係数を入力 ******
      e66(4,4)=gyz
      e66(5,5)=gxz
      e66(6,6)=gxy
!
!
      return
       end
