c 直交異方性材料の直方体構造物を、直方体要素の
c 線形／弾性の剛性行列で解くプログラム tyokum.f を
c 3種類の材料で解けるようにした。tyokuon.f で温度荷重を与える際、
c どうせ全節点に節点荷重を与えるのだから節点荷重ベクトルは、
c tyokuon.f の方で全節点の成分を定義したものを、そのまま
c 読み込むようにした。
c
c 剛性行列は、
c  Robert J. Melosh, Strucutural Analysis of Solids,
c  Jornal of the Structural Division, ASCE, 
c  Vol. 89, No. ST4, 1963, pp. 205-223
c にあるものを用いた。但し、印刷ミスと思われる以下の箇所について、
c  p.209: K11(1,1)の+4d'(1,1)を+4d'(1,1)に
c  p.209: K11(6,5)の+2d'(5,5)を+2d(5,5)に
c  p.209: K11(8,5)の-4d'(5,5)を-4d(5,5)に
c  p.209: K22(3,1)の+4d(4,4)を+d(4,4)に
c  p.210: K33(5,4)の-d'(6,6)を+d'(6,6)に
c  p.211: K31左(8,1)の2を-2に
c  p.211: k31左(8,4)の-1を-2に
c それぞれ修正した（対応箇所についてはプログラム中にも注意を書いた）
c
c 03/5/6版
c tyokus.f
c のサブルーチン solve をメインに入れて、配列サイズをもう少し
c 大きくとれるようにした。
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c プログラムを作った人：後藤文彦
c このプログラムは無保証です。
c このプログラムは改造／再配布して構いません。
c 但し、私的利用に留まらない場合は、もとのプログラムが
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/
c から取ってきたものと分かるように明示して下さい。
c 改造／再配布後も上記の条件を保持して下さい。
c バグの報告は、
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg
c にお願いします。
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
           program tyokuon3
           implicit real*8 (a-h, o-z)
       dimension bc(9400*3),f(9400*3),e66a(6,6),e66b(6,6),e66c(6,6),
     & d(9400*3),ngame(100),nfair(100),ijkab(9400),
     & kkaisi(9400),kryou(9400),bu(9400),bv(9400),bw(9400),
     & sm1(3,3,8,8),sm2(3,3,8,8),sm3(3,3,8,8),sm4(3,3,8,8),
     & sm5(3,3,8,8),sm6(3,3,8,8),n8(8),imx(9400*3),
     & su(3000000), sl(3000000),mksky(3000000),mkx(9400),
     & msl(9400*3),yokoza(100),tateza(100)
c 配列の最大サイズ：nmax
         nmax=3000000
c
c スカイラインの1次元配列のサイズは、同じ自由度数でも要素分割
c に依存するので、自由度数から一意には決められない。
c mage.f の要素分割のやり方だと、軸長（z軸）の分割数がxy断面
c の分割数に比べて多いとスパース部分が増える。
c 仮にスパース部分が皆無だったとすると
c およそ自由度×自由度／２のサイズになるが、こんなサイズは想定外。
c 取り敢えず、9400を最大節点数、その３倍を最大自由度数にしておく。
c 具体的な要素分割例としては、(縦+1)×(横+1)×(軸長+1)が、
c (6+1)*(13+1)*(100+1)=9191,(20+1)*(20+1)*(20+1)=9261 など。
c mage.f で幅（半解析だから幅の半分）×高さ×軸長が
c 1*(8～20)*(50～100)ぐらいの範囲のものを解く場合、有効数字３桁
c の収束解が得られる要素分割数は 6分割×10分割×100分割ぐらい。
c 幅（の半分）に対して桁高が10倍以上高い場合は、幅の分割数は
c 1 でもいいぐらいかも。
c
c 断面横方向の要素分割数：nx
         read(*,*) nx
c 断面縦方向の要素分割数：ny
         read(*,*) ny
c 軸方向の要素分割数：nz
         read(*,*) nz
c 要素数：nyou
            nyou=nx*ny*nz
c
c 要素番号 i の要素の分類番号
c yokoa*tateaの材料a：1
c yokoc*tateaの材料c：2
c yokob*tateaの材料b：3
c yokoa*tatebの材料b：4
c yokoc*tatebの材料b：5
c yokob*tatebの材料b：6
         do 13 i=1,nyou
          read(*,*) ijkab(i)
13       continue
c
c 材料aの要素の断面の横の長さ:yokoa, 縦の長さ:tatea
         read(*,*) yokoa
         read(*,*) tatea
c それらの半分
         yokoa1=yokoa/2.d0
         tatea1=tatea/2.d0
c
c 材料bの要素の断面の横の長さ:yokob, 縦の長さ:tateb
         read(*,*) yokob
         read(*,*) tateb
c それらの半分
         yokob1=yokob/2.d0
         tateb1=tateb/2.d0
c
c 材料cの要素の断面の横の長さ:yokoc, 縦の長さ:tatec
         read(*,*) yokoc
         read(*,*) tatec
c それらの半分
         yokoc1=yokoc/2.d0
         tatec1=tatec/2.d0
c
c 梁の軸長：ziku
         read(*,*) ziku
c 一要素の軸長の半分：ziku1
         ziku1=ziku/float(nz)/2.d0
c
c 節点数：nset
            nset=(nx+1)*(ny+1)*(nz+1)
c 自由度数：nziyuu
            nziyuu=nset*3
c
c ****** 境界条件 *******
       do 100 i=1, nziyuu
            bc(i)=1.d0
100    continue
           read(*,*) kkasyo
       do 110 i=1,kkasyo
            read(*,*) kkaisi(i)
            read(*,*) kryou(i)
            read(*,*) bu(i)
            read(*,*) bv(i)
            read(*,*) bw(i)
        do 110 j=kkaisi(i)-1, kryou(i)-1
            bc(3*j+1)=bu(i)
            bc(3*j+2)=bv(i)
            bc(3*j+3)=bw(i)
110    continue
c ****** 載荷条件 *******
       do 120 i=1, nziyuu
            f(i)=0.d0
120    continue
c 荷重ベクトル：f
        do 130 j=0, nset-1
          read(*,*) f(3*j+1)
          read(*,*) f(3*j+2)
          read(*,*) f(3*j+3)
130    continue
c
c 材料aの応力－ひずみ行列：e66a
c まずはe66aの3*3のところに「ひずみ－応力行列」を読み込む
          call zeromx(6,e66a)
         do i=1,3
         do j=1,3
            read(*,*) e66a(i,j)
         end do
         end do
           write(*,*) 'ea11=',e66a(1,1)
c で、ヤング率やポアソン比を書き写してから
            exxa=1.d0/e66a(1,1)
            eyya=1.d0/e66a(2,2)
            ezza=1.d0/e66a(3,3)
            poixya=-e66a(1,2)*exxa
            poixza=-e66a(1,3)*exxa
            poiyza=-e66a(2,3)*eyya
c 逆行列を取って、応力－ひずみ行列に変換する
          call gyaku(e66a)
c せん断弾性係数を読み込む
            read(*,*) e66a(4,4)
            read(*,*) e66a(5,5)
            read(*,*) e66a(6,6)
c 一応、せん断弾性係数も書き写しておく
            gxya=e66a(4,4)
            gxza=e66a(5,5)
            gyza=e66a(6,6)
c
c 材料bの応力－ひずみ行列：e66b
c まずはe66bの3*3のところに「ひずみ－応力行列」を読み込む
          call zeromx(6,e66b)
         do i=1,3
         do j=1,3
            read(*,*) e66b(i,j)
         end do
         end do
           write(*,*) 'eb11=',e66b(1,1)
c で、ヤング率やポアソン比を書き写してから
            exxb=1.d0/e66b(1,1)
            eyyb=1.d0/e66b(2,2)
            ezzb=1.d0/e66b(3,3)
            poixyb=-e66b(1,2)*exxb
            poixzb=-e66b(1,3)*exxb
            poiyzb=-e66b(2,3)*eyyb
c 逆行列を取って、応力－ひずみ行列に変換する
          call gyaku(e66b)
c せん断弾性係数を読み込む
            read(*,*) e66b(4,4)
            read(*,*) e66b(5,5)
            read(*,*) e66b(6,6)
c 一応、せん断弾性係数も書き写しておく
            gxyb=e66b(4,4)
            gxzb=e66b(5,5)
            gyzb=e66b(6,6)
c
c 材料cの応力－ひずみ行列：e66c
c まずはe66cの3*3のところに「ひずみ－応力行列」を読み込む
          call zeromx(6,e66c)
         do i=1,3
         do j=1,3
            read(*,*) e66c(i,j)
         end do
         end do
           write(*,*) 'ec11=',e66c(1,1)
c で、ヤング率やポアソン比を書き写してから
            exxc=1.d0/e66c(1,1)
            eyyc=1.d0/e66c(2,2)
            ezzc=1.d0/e66c(3,3)
            poixyc=-e66c(1,2)*exxc
            poixzc=-e66c(1,3)*exxc
            poiyzc=-e66c(2,3)*eyyc
c 逆行列を取って、応力－ひずみ行列に変換する
          call gyaku(e66c)
c せん断弾性係数を読み込む
            read(*,*) e66c(4,4)
            read(*,*) e66c(5,5)
            read(*,*) e66c(6,6)
c 一応、せん断弾性係数も書き写しておく
            gxyc=e66c(4,4)
            gxzc=e66c(5,5)
            gyzc=e66c(6,6)
c
c
c 荷重増分：dp
            read(*,*) dp
c
c 出力条件：
c 比較のために出力したい梁－柱理論値：nhikak
c nhikak=1: w 方向の伸びを出力
c
             read(*,*) nhikak
             write(*,*) 'nhikak=',nhikak
c 変位を画面出力したい節点の数：ngamen
             read(*,*) ngamen
c 変位を画面出力したい節点の節点番号：ngame(i)
         do 200 i=1,ngamen
             read(*,*) ngame(i)
200      continue
c 変位をファイル出力したい節点の数：nfairu
             read(*,*) nfairu
c 変位をファイル出力したい節点の節点番号：nfair(i)
         do 210 i=1,nfairu
             read(*,*) nfair(i)
210      continue
c
c 増分ステップ数：nstep
             read(*,*) nstep
c
c ***********データ入力終わり***********
             write(*,*) '材料aの応力-ひずみ行列：'
             write(*,212) ((e66a(i,j),j=1,3),i=1,3)
             write(*,*) 'Gxy,Gxz,Gyz='
             write(*,212) gxya,gxza,gyza
c
             write(*,*) '材料bの応力-ひずみ行列：'
             write(*,212) ((e66b(i,j),j=1,3),i=1,3)
             write(*,*) 'Gxy,Gxz,Gyz='
             write(*,212) gxyb,gxzb,gyzb
c
             write(*,*) '材料cの応力-ひずみ行列：'
             write(*,212) ((e66c(i,j),j=1,3),i=1,3)
             write(*,*) 'Gxy,Gxz,Gyz='
             write(*,212) gxyc,gxzc,gyzc
c
212         format(1p3d11.3)
c
c
c Melesh(1963)の剛性行列をsm1-sm6に書き込む
          call melesh(yokoa1,tatea1,ziku1,e66a,sm1)
          call melesh(yokoc1,tatea1,ziku1,e66c,sm2)
          call melesh(yokob1,tatea1,ziku1,e66b,sm3)
          call melesh(yokoa1,tateb1,ziku1,e66b,sm4)
          call melesh(yokoc1,tateb1,ziku1,e66b,sm5)
          call melesh(yokob1,tateb1,ziku1,e66b,sm6)
c
c まずは、スカイラインの概形をmkskyに書き込む
           mkx(1)=1
        do 980 i=1,nset-1
           mkx(i+1)=mkx(i)+i
980     continue
          lasto=mkx(nset)+nset-1
           write(*,*) ' 節点数(<9400)=',nset
           write(*,*) ' 自由度数(<28200)=',nziyuu
           write(*,*) ' スカイライン整形配列サイズ(<3000000)=',lasto
         if(lasto.gt.nmax) then
           write(*,*) ' スカイライン整形配列が大きすぎます'
           stop
         endif
c
          do 990 i=1,nset
          do 990 j=1,nset
             ij=mkx(j)+j-i
             mksky(ij)=0
990        continue
c
        do 1000 k=1,nz
        do 1000 j=1,ny
        do 1000 i=1,nx
         n8(8)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
         n8(4)=n8(8)+1
         n8(7)=n8(8)+(nx+1)
         n8(3)=n8(7)+1
         n8(5)=n8(8)+(nx+1)*(ny+1)
         n8(1)=n8(4)+(nx+1)*(ny+1)
         n8(6)=n8(7)+(nx+1)*(ny+1)
         n8(2)=n8(3)+(nx+1)*(ny+1)
c
        do 1000 ii=1,8
        do 1000 jj=1,8
           n8ii=n8(ii)
           n8jj=n8(jj)
          if(n8ii.gt.n8jj) goto 1000
           ij=mkx(n8jj)+n8jj-n8ii
           mksky(ij)=1
1000    continue
c
c スカイラインの対角項番号をimxに書き込む
        imx(1)=1
        do 1001 j=1,nset
           ifzumi=0
        do 1001 i=1,j
         ij=mkx(j)+j-i
         if(ifzumi.eq.1) goto 1004
          if(mksky(ij).eq.1) then
           ifzumi=1
           do 1002 k=1,3
            n3=3*(j-1)+k
            imx(n3+1)=imx(n3)+(j-i)*3+k
1002       continue
          endif
1004     continue
          if(ifzumi.eq.1) mksky(ij)=1
1001     continue
c
           imxyz=imx(nziyuu)+nziyuu-1
           write(*,*) ' スカイライン配列サイズ(<3000000)=  約',imxyz
          if(imxyz.ge.nmax) then
           write(*,*) 'たぶん配列サイズが大きすぎます'
           stop
          endif
c
c **************************************************************
             p=0.d0
        open(7,file='gnu.out')
        open(8,file='tanbu.out')
        open(9,file='ziku00.out')
        open(10,file='ziku10.out')
        open(11,file='ziku01.out')
        open(12,file='ziku11.out')
c 指定された荷重増分ステップまで以下の作業を繰り返す
        do 2000 i2000=1,nstep
             p=p+dp
c
c
c 全要素のsmを、それぞれの要素の8節点の節点番号ごとに
c 割り振って上三角スカイラインsu,下三角スカイラインslに書き込む
          do 220 i=1, imxyz
           su(i)=0.d0
           sl(i)=0.d0
220       continue
c
        do 1003 k=1,nz
        do 1003 j=1,ny
        do 1003 i=1,nx
         n8(8)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
         n8(4)=n8(8)+1
         n8(7)=n8(8)+(nx+1)
         n8(3)=n8(7)+1
         n8(5)=n8(8)+(nx+1)*(ny+1)
         n8(1)=n8(4)+(nx+1)*(ny+1)
         n8(6)=n8(7)+(nx+1)*(ny+1)
         n8(2)=n8(3)+(nx+1)*(ny+1)
c
        do 1003 ii=1,8
        do 1003 jj=1,8
           n8ii=n8(ii)
           n8jj=n8(jj)
        do 1003 l=1,3
        do 1003 m=1,3
           i24=(n8ii-1)*3 +l
           j24=(n8jj-1)*3 +m
        if(i24.le.j24) then
            ij=imx(j24)+j24-i24
            ijk=nx*ny*(k-1)+nx*(j-1)+i
            n1234=ijkab(ijk)
         goto (201,202,203,204,205,206), n1234
201       continue
            su(ij)=su(ij)+sm1(l,m,ii,jj)
          goto 207
202       continue
            su(ij)=su(ij)+sm2(l,m,ii,jj)
          goto 207
203       continue
            su(ij)=su(ij)+sm3(l,m,ii,jj)
          goto 207
204       continue
            su(ij)=su(ij)+sm4(l,m,ii,jj)
          goto 207
205       continue
            su(ij)=su(ij)+sm5(l,m,ii,jj)
          goto 207
206       continue
            su(ij)=su(ij)+sm6(l,m,ii,jj)
207      continue
        else
            ij=imx(i24)+i24-j24
            ijk=nx*ny*(k-1)+nx*(j-1)+i
            n1234=ijkab(ijk)
         goto (301,302,303,304,305,306), n1234
301       continue
            sl(ij)=sl(ij)+sm1(l,m,ii,jj)
          goto 307
302       continue
            sl(ij)=sl(ij)+sm2(l,m,ii,jj)
          goto 307
303       continue
            sl(ij)=sl(ij)+sm3(l,m,ii,jj)
          goto 307
304       continue
            sl(ij)=sl(ij)+sm4(l,m,ii,jj)
          goto 307
305       continue
            sl(ij)=sl(ij)+sm5(l,m,ii,jj)
          goto 307
306       continue
            sl(ij)=sl(ij)+sm6(l,m,ii,jj)
307      continue
        endif
1003    continue
c
c
         write(*,*) '境界条件'
c 境界条件を入れる
c su,slの中で定義されていないs(i,j)の空白部分を定義してはいけない
c 対角項以外
        do 1100 i=1,nset
        do 1100 j=i+1,nset
         mij=mkx(j)+j-i
         if(mksky(mij).eq.1) then
          do 1103 l=1,3
          do 1103 m=1,3
           i3=(i-1)*3+l
           j3=(j-1)*3+m
           ij=imx(j3)+j3-i3
           su(ij)=su(ij)*bc(i3)*bc(j3)
           sl(ij)=sl(ij)*bc(i3)*bc(j3)
1103      continue
         endif
1100    continue
c 対角項
        do 1101 i=1,nset
        do 1101 l=1,3
           i3=(i-1)*3+l
c 温度荷重などを全節点に与えてしまった時のため、
c 拘束節点の荷重ベクトルはここで0にする
           f(i3)=bc(i3)*f(i3)
        do 1101 m=l,3
           j3=(i-1)*3+m
           ij=imx(j3)+j3-i3
           su(ij)=su(ij)*bc(i3)*bc(j3)
           sl(ij)=sl(ij)*bc(i3)*bc(j3)
1101    continue
c
c
        do 1200 i=1,nziyuu
         if( bc(i).lt.1.d-3) then
            su(imx(i))=1.d0
            f(i)=0.d0
          end if
1200    continue
c
           do 2010 i=1,nziyuu
c 荷重項はslの対角項に入れる
            sl(imx(i))=f(i)*p
2010       continue
c
c
              write(*,*) '剛性方程式'
c 剛性方程式を解く
c tyokus.f
c (http://gthmhk.virtualave.net/programoj/tyoku/tyokus.html)
c では subroutine solve(stuh,stlh,imx,n01) というサブルーチン
c に置いていたものを、配列サイズ節約のため、メインに入れた。
c このもとのサブルーチンは、
c 岩熊先生(http://www.civil.tohoku.ac.jp/~bear/）が、
c 東京大学の大計のライブラリーかどこかの「ガウスの消去法」の
c サブルーチンを2次元骨組をスカイライン法＋弧長法で解くために
c 書き換えたものを私が3次元梁用に書き換えたもの
c (http://gthmhk.virtualave.net/programoj/hari/hari.html)
c を弧長法なし（＋固有値の計算なし）に書き換えた。
c
c
c 方程式は掃き出し法によりここで解かれる
c
c
         n01=nziyuu
c
        do 101 i=1,n01
        msl(i)=1+i+imx(i)-imx(i+1)
101     continue
c
c 前進消去
c
      do 10 j=2,n01
      ist=msl(j)+1
      ied=j-1
      if( ist.gt.ied ) go to 20
      do 30 i=ist,ied
      kst=msl(j)
      if( msl(i).gt.kst ) kst=msl(i)
      ked=i-1
      ij=imx(j)+j-i
      do 30 k=kst,ked
      kj=imx(j)+j-k
      ki=imx(i)+i-k
      su(ij)=su(ij)-su(kj)*sl(ki)
      sl(ij)=sl(ij)-sl(kj)*su(ki)
   30 continue
   20 continue
      do 40 i=msl(j),ied
      ij=imx(j)+j-i
      su(ij)=su(ij)/su(imx(i))
      sl(ij)=sl(ij)/su(imx(i))
   40 continue
      do 10 k=msl(j),ied
      kj=imx(j)+j-k
      su(imx(j))=su(imx(j))-su(kj)*sl(kj)*su(imx(k))
      sl(imx(j))=sl(imx(j))-sl(kj)*sl(imx(k))
   10 continue
c
c 後進消去
c
      do 60 j=1,n01
      sl(imx(j))=sl(imx(j))/su(imx(j))
   60 continue
      do 70 j=2,n01
      i=n01+2-j
      ked=i-1
      do 70 k=msl(i),ked
      ii=imx(i)+i-k
      sl(imx(k))=sl(imx(k))-su(ii)*sl(imx(i))
   70 continue
c     return
c      end
c
c
c
c
        do i=1,nziyuu
         d(i)=sl(imx(i))
        end do
c
c
c ************ 出力 ***************
c
       do 3000 j=1,ngamen
        nga=ngame(j)*3-3
         write(*,*) '--------------------------------------------------'
         write(*,*) '節点番号：',ngame(j)
c
         aa=2.d0*yokoa*tatea*2.d0
         ea=ezza*aa+ezzb*(a-aa)
         yoko2=(float(nx)-1.d0)*yokob
         yoko=yokoa+yoko2
         tate2=(float(ny)-4.d0)*tateb
         tate=4.d0*tatea+tate2
         a=yoko*tate
         xi=yoko*tate**3/12.d0
         xia=(yokoa*tate**3-yokoa*tate2**3)/12.d0
         xib=xi-xia
         ezzxi=ezza*xia+ezzb*xib
         yi=tate*yoko**3/12.d0
c
          goto 599
          goto (511,512,513), nhikak
511    continue
c        write(*,125) ' f=',p,'   PL/EA=',p*ziku/ea
       goto 599
512    continue
        tawa=p*ziku**3/3.d0/ezzxi
        pk=1.2d0
        tawap=tawa+p*ziku/gyzb/a*pk
       write(*,115) ' 荷重 f=',p
       write(*,115) ' たわみ v(せん断無視の梁理論)=',tawa
       write(*,115) ' たわみ v(νから  せん断考慮)=',tawap
114    format(a, f10.3)
       goto 599
112    format(1p2d13.5)
513    write(*,125) ' f=',p,'   PL^3/3EI=',p*ziku**3/3.d0/ezz/yi
       goto 599
         write(*,115) ' f=',p
599     continue
         write(*,135) ' u,v,w=',d(nga+1),d(nga+2),d(nga+3)
3000   continue
c
c
c 自由端部の全節点変位をgnuplotの３次元プロット用に出力
        yokoza(1)=0.d0
         j=1
         k=nz
        do i=1,nx
         ijk=nx*ny*(k-1)+nx*(j-1)+i
         n1234=ijkab(ijk)
         if(n1234.eq.1) then
          yokoza(i+1)=yokoa
         else
          if(n1234.eq.2) then
           yokoza(i+1)=yokoc
          else
           yokoza(i+1)=yokob
          end if
         end if
        end do
c
        do i=2,nx+1
         yokoza(i)=yokoza(i)+yokoza(i-1)
        end do
c
        tateza(1)=0.d0
         i=1
         k=nz
        do j=1,ny
         ijk=nx*ny*(k-1)+nx*(j-1)+i
         n1234=ijkab(ijk)
         if(n1234.eq.1) then
          tateza(j+1)=tatea
         else
          tateza(j+1)=tateb
         end if
        end do
c
        do j=2,ny+1
         tateza(j)=tateza(j)+tateza(j-1)
        end do
c
        ntan0=(nx+1)*(ny+1)*nz
       do j=1,ny+1
        i1=ntan0+(nx+1)*(j-1)
        do i=1,nx+1
         i2=i1+i
c 変位に倍率をかける時は、下の「*1.d0」を適当に変える
         dx=d(i2*3-2)* 1.d0
         dy=d(i2*3-1)* 1.d0
         dz=d(i2*3)  * 1.d0
        write(7,185) dx+yokoza(i),dy+tateza(j),dz
        end do
        write(7,*)
       end do
c
c 自由端部の全節点変位の出力
        nfair1=nfair(1)
        write(8,*) '-------------------------------------------------'
        write(8,*) 'u'
        do i=1,ny+1
         j1=nfair(1)+(nx+1)*(i-1)
         write(8,189) (d(j*3-2),j=j1,j1+nx)
        end do
        write(8,*) '-------------------------------------------------'
        write(8,*) 'v'
        do i=1,ny+1
         j1=nfair(1)+(nx+1)*(i-1)
         write(8,189) (d(j*3-1),j=j1,j1+nx)
        end do
        write(8,*) '-------------------------------------------------'
        write(8,*) 'w'
        do i=1,ny+1
         j1=nfair(1)+(nx+1)*(i-1)
         write(8,189) (d(j*3),j=j1,j1+nx)
        end do
        write(8,*) '-------------------------------------------------'
c
c 梁の4角の4辺の全節点変位をそれぞれ軸方向に出力
        do k=1,nz+1
          za=ziku1*2.d0*float(k-1)
           n00=(nx+1)*(ny+1)*(k-1)+1
           n10=(nx+1)*(ny+1)*(k-1)+nx+1
           n01=(nx+1)*(ny+1)*(k-1)+(nx+1)*ny+1
           n11=(nx+1)*(ny+1)*k
           write( 9,195) za,d(n00*3-2),d(n00*3-1),d(n00*3)
           write(10,195) za,d(n10*3-2),d(n10*3-1),d(n10*3)
           write(11,195) za,d(n01*3-2),d(n01*3-1),d(n01*3)
           write(12,195) za,d(n11*3-2),d(n11*3-1),d(n11*3)
      end do
c
2000   continue
        close(7)
        close(8)
        close(9)
        close(10)
        close(11)
        close(12)
115      format(a,1pd13.5)
125      format(a,1pd12.5,a,1pd12.5)
126      format(a,1pd12.5,a,1pd12.5,a,1pd12.5)
135      format(a,1p3d13.5)
145      format(1p4d13.5)
185      format(1p3d13.5)
189      format(1p9d11.3)
175      format(1p7d13.5)
195      format(1p4d13.5)
c
c ouryoku.fで各要素の応力を求める時のために節点変位dを
c formatなしでfort.12に書き出しておく。
        write(12) d
        rewind(12)
c
           end
c
c
c **********************************************************
c
c ひずみ－応力行列の逆数を取って、
c 応力－ひずみ行列を求める。
c
            subroutine gyaku(a)
             implicit real*8(a-h, o-z)
             dimension a(6,6), b(3,3)
              det=a(1,1)*a(2,2)*a(3,3)+a(2,1)*a(3,2)*a(1,3)
     &           +a(3,1)*a(1,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)
     &           -a(3,1)*a(2,2)*a(1,3)-a(1,1)*a(3,2)*a(2,3)
              b(1,1)=(-1.d0)**(1+1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))
              b(1,2)=(-1.d0)**(1+2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))
              b(1,3)=(-1.d0)**(1+3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2))
              b(2,1)=(-1.d0)**(2+1)*(a(1,2)*a(3,3)-a(3,2)*a(1,3))
              b(2,2)=(-1.d0)**(2+2)*(a(1,1)*a(3,3)-a(3,1)*a(1,3))
              b(2,3)=(-1.d0)**(2+3)*(a(1,1)*a(3,2)-a(3,1)*a(1,2))
              b(3,1)=(-1.d0)**(3+1)*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
              b(3,2)=(-1.d0)**(3+2)*(a(1,1)*a(2,3)-a(2,1)*a(1,3))
              b(3,3)=(-1.d0)**(3+3)*(a(1,1)*a(2,2)-a(2,1)*a(1,2))
             do 10 i=1,3
             do 10 j=1,3
              a(i,j)=b(j,i)/det
10           continue
             return
             end
c
c
c 剛性行列
        subroutine melesh(a,b,c,d,sm)
        implicit real*8(a-h, o-z)
         dimension d(6,6),sm(3,3,8,8),
     & s21a(8,8),s21b(8,8),s31a(8,8),s31b(8,8),s32a(8,8),s32b(8,8)
c sm に Melesh(1963) の剛性行列そのままを与えて
c メインプログラムで節点順に並び替えた剛性行列を s に入れる。
          d11=d(1,1)*b**2*c**2
          d22=d(2,2)*a**2*c**2
          d33=d(3,3)*a**2*b**2
          d44=d(4,4)*b**2*c**2
          d55=d(5,5)*a**2*b**2
          d66=d(6,6)*a**2*b**2
          d44d=d(4,4)*a**2*c**2
          d55d=d(5,5)*b**2*c**2
          d66d=d(6,6)*a**2*c**2
          d44dd=c*d(4,4)
          d55dd=b*d(5,5)
          d66dd=a*d(6,6)
          d12=c*d(1,2)
          d13=b*d(1,3)
          d23=a*d(2,3)
c
          sm(1,1,1,1)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,2,1)= 2.d0*d11 -4.d0*d44d+2.d0*d55
          sm(1,1,3,1)=      d11 -2.d0*d44d-2.d0*d55
          sm(1,1,4,1)= 2.d0*d11 +2.d0*d44d-4.d0*d55
          sm(1,1,5,1)=-4.d0*d11 +2.d0*d44d+2.d0*d55
          sm(1,1,6,1)=-2.d0*d11 -2.d0*d44d     +d55
          sm(1,1,7,1)=     -d11      -d44d     -d55
          sm(1,1,8,1)=-2.d0*d11      +d44d-2.d0*d55
c Melesh(63)ではsm(1,1,1,1)のd44dがd44になっている。
c
          sm(1,1,2,2)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,3,2)= 2.d0*d11 +2.d0*d44d-4.d0*d55
          sm(1,1,4,2)=      d11 -2.d0*d44d-2.d0*d55
          sm(1,1,5,2)=-2.d0*d11 -2.d0*d44d     +d55
          sm(1,1,6,2)=-4.d0*d11 +2.d0*d44d+2.d0*d55
          sm(1,1,7,2)=-2.d0*d11      +d44d-2.d0*d55
          sm(1,1,8,2)=     -d11      -d44d     -d55
c
          sm(1,1,3,3)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,4,3)= 2.d0*d11 -4.d0*d44d+2.d0*d55
          sm(1,1,5,3)=     -d11      -d44d     -d55
          sm(1,1,6,3)=-2.d0*d11      +d44d-2.d0*d55
          sm(1,1,7,3)=-4.d0*d11 +2.d0*d44d+2.d0*d55
          sm(1,1,8,3)=-2.d0*d11 -2.d0*d44d     +d55
c
          sm(1,1,4,4)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,5,4)=-2.d0*d11      +d44d-2.d0*d55
          sm(1,1,6,4)=     -d11      -d44d     -d55
          sm(1,1,7,4)=-2.d0*d11 -2.d0*d44d     +d55
          sm(1,1,8,4)=-4.d0*d11 +2.d0*d44d+2.d0*d55
c
          sm(1,1,5,5)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,6,5)= 2.d0*d11 -4.d0*d44d+2.d0*d55
          sm(1,1,7,5)=      d11 -2.d0*d44d-2.d0*d55
          sm(1,1,8,5)= 2.d0*d11 +2.d0*d44d-4.d0*d55
c Melesh(1963)ではsm(1,1,6,5)とsm(1,1,8,5)のd55がd55dになっている
c
          sm(1,1,6,6)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,7,6)= 2.d0*d11 +2.d0*d44d-4.d0*d55
          sm(1,1,8,6)=      d11 -2.d0*d44d-2.d0*d55
c
          sm(1,1,7,7)= 4.d0*d11 +4.d0*d44d+4.d0*d55
          sm(1,1,8,7)= 2.d0*d11 -4.d0*d44d+2.d0*d55
c
          sm(1,1,8,8)= 4.d0*d11 +4.d0*d44d+4.d0*d55
c
c
          sm(2,2,1,1)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,2,1)=-4.d0*d22 +2.d0*d44 +2.d0*d66
          sm(2,2,3,1)=-2.d0*d22      +d44 -2.d0*d66
          sm(2,2,4,1)= 2.d0*d22 +2.d0*d44 -4.d0*d66
          sm(2,2,5,1)= 2.d0*d22 -4.d0*d44 +2.d0*d66
          sm(2,2,6,1)=-2.d0*d22 -2.d0*d44      +d66
          sm(2,2,7,1)=     -d22      -d44      -d66
          sm(2,2,8,1)=      d22 -2.d0*d44 -2.d0*d66
c Melesh(63)では、sm(2,2,3,1)のd44が4.d0*d44になっている
c         sm(2,2,3,1)=-2.d0*d22 +4.d0*d44 -2.d0*d66
c
          sm(2,2,2,2)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,3,2)= 2.d0*d22 +2.d0*d44 -4.d0*d66
          sm(2,2,4,2)=-2.d0*d22      +d44 -2.d0*d66
          sm(2,2,5,2)=-2.d0*d22 -2.d0*d44      +d66
          sm(2,2,6,2)= 2.d0*d22 -4.d0*d44 +2.d0*d66
          sm(2,2,7,2)=      d22 -2.d0*d44 -2.d0*d66
          sm(2,2,8,2)=     -d22      -d44      -d66
c
          sm(2,2,3,3)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,4,3)=-4.d0*d22 +2.d0*d44 +2.d0*d66
          sm(2,2,5,3)=     -d22      -d44      -d66
          sm(2,2,6,3)=      d22 -2.d0*d44 -2.d0*d66
          sm(2,2,7,3)= 2.d0*d22 -4.d0*d44 +2.d0*d66
          sm(2,2,8,3)=-2.d0*d22 -2.d0*d44      +d66
c
          sm(2,2,4,4)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,5,4)=      d22 -2.d0*d44 -2.d0*d66
          sm(2,2,6,4)=     -d22      -d44      -d66
          sm(2,2,7,4)=-2.d0*d22 -2.d0*d44      +d66
          sm(2,2,8,4)= 2.d0*d22 -4.d0*d44 +2.d0*d66
c
          sm(2,2,5,5)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,6,5)=-4.d0*d22 +2.d0*d44 +2.d0*d66
          sm(2,2,7,5)=-2.d0*d22      +d44 -2.d0*d66
          sm(2,2,8,5)= 2.d0*d22 +2.d0*d44 -4.d0*d66
c
          sm(2,2,6,6)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,7,6)= 2.d0*d22 +2.d0*d44 -4.d0*d66
          sm(2,2,8,6)=-2.d0*d22      +d44 -2.d0*d66
c
          sm(2,2,7,7)= 4.d0*d22 +4.d0*d44 +4.d0*d66
          sm(2,2,8,7)=-4.d0*d22 +2.d0*d44 +2.d0*d66
c
          sm(2,2,8,8)= 4.d0*d22 +4.d0*d44 +4.d0*d66
c
c
          sm(3,3,1,1)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,2,1)= 2.d0*d33 +2.d0*d55d -4.d0*d66d
          sm(3,3,3,1)=-2.d0*d33      +d55d -2.d0*d66d
          sm(3,3,4,1)=-4.d0*d33 +2.d0*d55d +2.d0*d66d
          sm(3,3,5,1)= 2.d0*d33 -4.d0*d55d +2.d0*d66d
          sm(3,3,6,1)=      d33 -2.d0*d55d -2.d0*d66d
          sm(3,3,7,1)=     -d33      -d55d      -d66d
          sm(3,3,8,1)=-2.d0*d33 -2.d0*d55d      +d66d
c
          sm(3,3,2,2)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,3,2)=-4.d0*d33 +2.d0*d55d +2.d0*d66d
          sm(3,3,4,2)=-2.d0*d33      +d55d -2.d0*d66d
          sm(3,3,5,2)=      d33 -2.d0*d55d -2.d0*d66d
          sm(3,3,6,2)= 2.d0*d33 -4.d0*d55d +2.d0*d66d
          sm(3,3,7,2)=-2.d0*d33 -2.d0*d55d      +d66d
          sm(3,3,8,2)=     -d33      -d55d      -d66d
c
          sm(3,3,3,3)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,4,3)= 2.d0*d33 +2.d0*d55d -4.d0*d66d
          sm(3,3,5,3)=     -d33      -d55d      -d66d
          sm(3,3,6,3)=-2.d0*d33 -2.d0*d55d      +d66d
          sm(3,3,7,3)= 2.d0*d33 -4.d0*d55d +2.d0*d66d
          sm(3,3,8,3)=      d33 -2.d0*d55d -2.d0*d66d
c
          sm(3,3,4,4)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
c          sm(3,3,5,4)=-2.d0*d33 -2.d0*d55d      -d66d
          sm(3,3,5,4)=-2.d0*d33 -2.d0*d55d      +d66d
          sm(3,3,6,4)=     -d33      -d55d      -d66d
          sm(3,3,7,4)=      d33 -2.d0*d55d -2.d0*d66d
          sm(3,3,8,4)= 2.d0*d33 -4.d0*d55d +2.d0*d66d
c Melesh(63)ではsm(3,3,5,4)の+d66dは-d66dになっているが、
c 答えが出ないので+d66dに変えた。
c
          sm(3,3,5,5)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,6,5)= 2.d0*d33 +2.d0*d55d -4.d0*d66d
          sm(3,3,7,5)=-2.d0*d33      +d55d -2.d0*d66d
          sm(3,3,8,5)=-4.d0*d33 +2.d0*d55d +2.d0*d66d
c
          sm(3,3,6,6)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,7,6)=-4.d0*d33 +2.d0*d55d +2.d0*d66d
          sm(3,3,8,6)=-2.d0*d33      +d55d -2.d0*d66d
c
          sm(3,3,7,7)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
          sm(3,3,8,7)= 2.d0*d33 +2.d0*d55d -4.d0*d66d
c
          sm(3,3,8,8)= 4.d0*d33 +4.d0*d55d +4.d0*d66d
c
c
          do 10 n=1,3
          do 10 i=1,8
           do 15 j=1,i
            sm(n,n,i,j)=sm(n,n,i,j)/18.d0/a/b/c
15         continue
          do 10 j=1,i-1
            sm(n,n,j,i)=sm(n,n,i,j)
10        continue
c
c
           s21a(1,1)=-2.d0
           s21a(2,1)= 2.d0
           s21a(3,1)= 1.d0
           s21a(4,1)=-1.d0
           s21a(5,1)=-2.d0
           s21a(6,1)= 2.d0
           s21a(7,1)= 1.d0
           s21a(8,1)=-1.d0
c
           s21a(1,3)=-1.d0
           s21a(2,3)= 1.d0
           s21a(3,3)= 2.d0
           s21a(4,3)=-2.d0
           s21a(5,3)=-1.d0
           s21a(6,3)= 1.d0
           s21a(7,3)= 2.d0
           s21a(8,3)=-2.d0
c
           s21b(1,1)=-2.d0
           s21b(2,1)=-2.d0
           s21b(3,1)=-1.d0
           s21b(4,1)=-1.d0
           s21b(5,1)= 2.d0
           s21b(6,1)= 2.d0
           s21b(7,1)= 1.d0
           s21b(8,1)= 1.d0
c
           s21b(1,3)= 1.d0
           s21b(2,3)= 1.d0
           s21b(3,3)= 2.d0
           s21b(4,3)= 2.d0
           s21b(5,3)=-1.d0
           s21b(6,3)=-1.d0
           s21b(7,3)=-2.d0
           s21b(8,3)=-2.d0
c
           s31a(1,1)= 2.d0
           s31a(2,1)= 1.d0
           s31a(3,1)=-1.d0
           s31a(4,1)=-2.d0
           s31a(5,1)= 2.d0
           s31a(6,1)= 1.d0
           s31a(7,1)=-1.d0
           s31a(8,1)=-2.d0
c Melesh(63)ではs31a(8,1)は2になっているが、i=1,8まで
c 足しても0にならないので -2 の間違いだと思う。
c ということは、Melesh(63)でs31a(8,4)が -1 となっているのも
c -2 の間違いだと思う。
           s31a(1,2)= 1.d0
           s31a(2,2)= 2.d0
           s31a(3,2)=-2.d0
           s31a(4,2)=-1.d0
           s31a(5,2)= 1.d0
           s31a(6,2)= 2.d0
           s31a(7,2)=-2.d0
           s31a(8,2)=-1.d0
c
           s31b(1,1)= 2.d0
           s31b(2,1)= 1.d0
           s31b(3,1)= 1.d0
           s31b(4,1)= 2.d0
           s31b(5,1)=-2.d0
           s31b(6,1)=-1.d0
           s31b(7,1)=-1.d0
           s31b(8,1)=-2.d0
c
           s31b(1,2)= 1.d0
           s31b(2,2)= 2.d0
           s31b(3,2)= 2.d0
           s31b(4,2)= 1.d0
           s31b(5,2)=-1.d0
           s31b(6,2)=-2.d0
           s31b(7,2)=-2.d0
           s31b(8,2)=-1.d0
c
           s32a(1,1)=-2.d0
           s32a(2,1)=-2.d0
           s32a(3,1)= 2.d0
           s32a(4,1)= 2.d0
           s32a(5,1)=-1.d0
           s32a(6,1)=-1.d0
           s32a(7,1)= 1.d0
           s32a(8,1)= 1.d0
c
           s32a(1,5)=-1.d0
           s32a(2,5)=-1.d0
           s32a(3,5)= 1.d0
           s32a(4,5)= 1.d0
           s32a(5,5)=-2.d0
           s32a(6,5)=-2.d0
           s32a(7,5)= 2.d0
           s32a(8,5)= 2.d0
c
           s32b(1,1)=-2.d0
           s32b(2,1)= 2.d0
           s32b(3,1)= 2.d0
           s32b(4,1)=-2.d0
           s32b(5,1)=-1.d0
           s32b(6,1)= 1.d0
           s32b(7,1)= 1.d0
           s32b(8,1)=-1.d0
c
           s32b(1,5)=-1.d0
           s32b(2,5)= 1.d0
           s32b(3,5)= 1.d0
           s32b(4,5)=-1.d0
           s32b(5,5)=-2.d0
           s32b(6,5)= 2.d0
           s32b(7,5)= 2.d0
           s32b(8,5)=-2.d0
c
           do 20 i=1,8
            s21a(i,2)= s21a(i,1)
            s21a(i,5)=-s21a(i,1)
            s21a(i,6)=-s21a(i,1)
            s21a(i,4)= s21a(i,3)
            s21a(i,7)=-s21a(i,3)
            s21a(i,8)=-s21a(i,3)
c
            s21b(i,2)=-s21b(i,1)
            s21b(i,5)= s21b(i,1)
            s21b(i,6)=-s21b(i,1)
            s21b(i,4)=-s21b(i,3)
            s21b(i,7)= s21b(i,3)
            s21b(i,8)=-s21b(i,3)
c
            s31a(i,4)= s31a(i,1)
            s31a(i,5)=-s31a(i,1)
            s31a(i,8)=-s31a(i,1)
            s31a(i,3)= s31a(i,2)
            s31a(i,6)=-s31a(i,2)
            s31a(i,7)=-s31a(i,2)
c
            s31b(i,4)=-s31b(i,1)
            s31b(i,5)= s31b(i,1)
            s31b(i,8)=-s31b(i,1)
            s31b(i,3)=-s31b(i,2)
            s31b(i,6)= s31b(i,2)
            s31b(i,7)=-s31b(i,2)
c
            s32a(i,2)=-s32a(i,1)
            s32a(i,3)=-s32a(i,1)
            s32a(i,4)= s32a(i,1)
            s32a(i,6)=-s32a(i,5)
            s32a(i,7)=-s32a(i,5)
            s32a(i,8)= s32a(i,5)
c
            s32b(i,2)= s32b(i,1)
            s32b(i,3)=-s32b(i,1)
            s32b(i,4)=-s32b(i,1)
            s32b(i,6)= s32b(i,5)
            s32b(i,7)=-s32b(i,5)
            s32b(i,8)=-s32b(i,5)
20         continue
c Melesh(63)ではs31a(8,4)=-1.d0となっているが、これだと
c (8,j),j=1,8まで総和したときに0にならない。
c
       do 30 i=1,8
       do 30 j=1,8
           sm(2,1,i,j)=( d12*s21a(i,j) + d44dd*s21b(i,j) )/12.d0
           sm(3,1,i,j)=( d13*s31a(i,j) + d55dd*s31b(i,j) )/12.d0
           sm(3,2,i,j)=( d23*s32a(i,j) + d66dd*s32b(i,j) )/12.d0
30     continue
       do 40 i=1,8
       do 40 j=1,8
           sm(1,2,i,j)=sm(2,1,j,i)
           sm(1,3,i,j)=sm(3,1,j,i)
           sm(2,3,i,j)=sm(3,2,j,i)
40     continue
c
           return
           end
c
c
        subroutine zeromx(n,a)
        implicit real*8(a-h,o-z)
         dimension a(n,n)
          do 10 i=1,n
          do 10 j=1,n
             a(i,j)=0.d0
10        continue
         return
        end
c
