c tyokuon.fで出力される節点変位にひずみ-変位行列と応力-ひずみ行列を
c かけて各要素の応力を求める。
c 手順は、
c hippou を実行してhippon.d を作る。
c このとき、ひずみ-変位行列や応力-ひずみ行列などがfort.11に出力される。
c tyokuon<hippon.d を実行すると変位dがfort.12に出力される。
c ouryokuを実行すると、fort.11とfort.12を読み込んで各要素の応力を求める。
c 03/3/17版
c
        program ouryoku
       implicit real*8 (a-h, o-z)
       dimension ijkab(9400),e66a(6,6),e66b(6,6),
     & d(9400*3),d24(24),n8(8),sgm(9400,6),sgm6(6),bd6(6),
     & bt1(24,6),bt2(24,6),bt3(24,6),bt4(24,6)
c
       read(11) nx,ny,nz,ijkab,e66a,e66b,bt1,bt2,bt3,bt4
       rewind(11)
       read(12) d
       rewind(12)
        write(*,*) 'nx,ny,nz=',nx,ny,nz
c 全節点数：nset
          nset=(nx+1)*(ny+1)*(nz+1)
c 全自由度数：nziyuu
c
c 要素ごとに8節点の節点番号を求め、その番号に対応する変位を
c d24に読み込む。
c
        do 1001 k=1,nz
        do 1001 j=1,ny
        do 1001 i=1,nx
         n8(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
         n8(2)=n8(1)+1
         n8(3)=n8(1)+(nx+1)
         n8(4)=n8(3)+1
         n8(5)=n8(1)+(nx+1)*(ny+1)
         n8(6)=n8(2)+(nx+1)*(ny+1)
         n8(7)=n8(3)+(nx+1)*(ny+1)
         n8(8)=n8(4)+(nx+1)*(ny+1)
c
        ijk=nx*ny*(k-1)+nx*(j-1)+i
        nab=ijkab(ijk)
c
        do  m=1,3
        do  jj=1,8
           j24=(n8(jj)-1)*3
           j83=(jj-1)*3+m
            d24(j83)=d(j24+m)
         end do
         end do
c
c 以下の分類番号に応じたひずみ-変位行列と応力-ひずみ行列を
c 要素ごとの節点変位にかける
c yokoa*tateaの鋼：nab=1
c yokob*tateaの木：nab=2
c yokoa*tatebの木：nab=3
c yokob*tatebの木：nab=4
       goto (401,402,403,404), nab
c
401    continue
        call bd(bt1,d24,bd6)
        call ebd(e66a,bd6,sgm6)
       goto 444
c
402    continue
        call bd(bt2,d24,bd6)
        call ebd(e66b,bd6,sgm6)
       goto 444
c
403    continue
        call bd(bt3,d24,bd6)
        call ebd(e66b,bd6,sgm6)
       goto 444
c
404    continue
        call bd(bt4,d24,bd6)
        call ebd(e66b,bd6,sgm6)
444    continue
c
c 要素応力を要素番号ごとにsgmに書き込む。
       do i6=1,6
        sgm(ijk,i6)=sgm6(i6)
       end do
c
1001    continue
c
c *********** 出力例 ******************
c 例えば固定端と自由端に面した全要素の応力を書き出す場合：
           k=1
        do j=1,ny
        do i=1,nx
         write(*,*) 'i,j,k=',i,j,k
         ijk=nx*ny*(k-1)+nx*(j-1)+i
         write(*,200) (sgm(ijk,n6),n6=1,6)
        end do
        end do
       if(nz.ne.1) then
           k=nz
        do j=1,ny
        do i=1,nx
         write(*,*) 'i,j,k=',i,j,k
         ijk=nx*ny*(k-1)+nx*(j-1)+i
         write(*,200) (sgm(ijk,n6),n6=1,6)
        end do
        end do
       end if
200    format(1p6d13.5)
c
       end
c
c
       subroutine bd(bt,d24,bd6)
       implicit real*8 (a-h, o-z)
       dimension bt(24,6),d24(24),bd6(6)
        do i=1,6
          bd6(i)=0.d0
          do j=1,24
            bd6(i)=bd6(i)+bt(j,i)*d24(j)
          end do
        end do
       return
       end
c
       subroutine ebd(e66,bd6,sgm6)
       implicit real*8 (a-h, o-z)
       dimension e66(6,6),bd6(6),sgm6(6)
        do i=1,6
          sgm6(i)=0.d0
          do j=1,6
            sgm6(i)=sgm6(i)+e66(i,j)*bd6(j)
          end do
        end do
       return
       end
c
