c 直交異方性材料の直方体構造物を、直方六面体（ソリッド）要素の
c 線形／弾性の剛性行列で解くプログラム。剛性行列は、
c Robert J. Melesh, 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(6,5)の+2d'(5,5)を+2d(5,5)に
c p.209: K11(6,8)の-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 02/5/24版
c プログラムを作った人：後藤文彦
c 改造などは自由ですが、もとのプログラムが、
c http://gthmhk.virtualave.net/programoj/tyoku/tyoku.f
c であることをプログラム中に明記して戴けると嬉しいです。
c バグの報告は、
c http://gthmhk.virtualave.net/denbin.html
c にお願いします。
c
           program tyoku
           implicit real*8 (a-h, o-z)
       dimension bc(2700*3),f(2700*3),e66(6,6),
     & d(2700*3),ngame(100),nfair(100),
     & kkaisi(2700),kryou(2700),bu(2700),bv(2700),bw(2700),
     & lkaisi(2700),lryou(2700),fx(2700),fy(2700),fz(2700),
     & sm(3,3,8,8),s(2700*3,2700*3),n8(8),
     & st(2700*3,2700*3),df(2700*3)
c
c 入力可能な自由度の最大値：nmax
         nmax=2700*3
c 私の環境（cygwin の g77）だと配列の最大サイズは a(67108863)
c つまり 2**26=67108864 未満。
c tyoku.f では自由度×自由度の行列を解くので、
c (2700*3)*(2700*3)=65610000 を最大サイズにしておく。
c 具体的な要素分割例としては、縦×横×軸長が、
c 6*6*75=2700, 8*8*42=2688<2700, 13*13*13=2199<2700 など。
c
c 梁の断面の横の長さ:yoko, 縦の長さ:tate
         read(*,*) yoko
         read(*,*) tate
c 梁の軸長：ziku
         read(*,*) ziku
c
         a=yoko*tate
         xi=yoko*tate**3/12.d0
         yi=tate*yoko**3/12.d0
c 断面横方向の要素分割数：nx
         read(*,*) nx
c 一要素の横の長さの半分:yoko1
         yoko1=yoko/float(nx)/2.d0
c 断面縦方向の要素分割数：ny
         read(*,*) ny
c 一要素の縦の長さの半分：tate1
         tate1=tate/float(ny)/2.d0
c 軸方向の要素分割数：nz
         read(*,*) nz
c 一要素の軸長の半分：ziku1
         ziku1=ziku/float(nz)/2.d0
c
c 要素数：nyou
            nyou=nx*ny*nz
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
            read(*,*) lkasyo
       do 130 i=1,lkasyo
            read(*,*) lkaisi(i)
            read(*,*) lryou(i)
            read(*,*) fx(i)
            read(*,*) fy(i)
            read(*,*) fz(i)
        do 130 j=lkaisi(i)-1, lryou(i)-1
            f(3*j+1)=fx(i)
            f(3*j+2)=fy(i)
            f(3*j+3)=fz(i)
130    continue
c 応力－ひずみ行列：e66
c 最初は、e66の3*3のところに、ひずみ－応力行列を読み込む
          call zeromx(6,e66)
         do 140 i=1,3
         do 140 j=1,3
            read(*,*) e66(i,j)
140      continue
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 逆行列を取って、応力－ひずみ行列に変換する
          call gyaku(e66)
c せん断弾性係数を読み込む
            read(*,*) e66(4,4)
            read(*,*) e66(5,5)
            read(*,*) e66(6,6)
c 一応、せん断弾性係数も書き写しておく
            gxy=e66(4,4)
            gxz=e66(5,5)
            gyz=e66(6,6)
c 荷重増分：dp
            read(*,*) dp
c
c 出力条件：
c 比較のために出力したい梁－柱理論値：nhikak
c nhikak=1: w 方向の伸びを出力
c
             read(*,*) 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 ***********データ入力終わり***********
c
c
c Melesh(1963)の剛性行列をsmに書き込む
          call melesh(yoko1,tate1,ziku1,e66,sm)
c
          call zeromx(nziyuu,s)
c
c 全要素のsmを、それぞれの要素の8節点の節点番号ごとに
C 割り振ってsに書き込む
        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
c
        do 1000 l=1,3
        do 1000 m=1,3
        do 1000 ii=1,8
        do 1000 jj=1,8
           i24=(n8(ii)-1)*3
           j24=(n8(jj)-1)*3
            s(i24+l,j24+m)=s(i24+l,j24+m)+sm(l,m,ii,jj)
1000    continue
c
c
c 境界条件を入れる
        do 1100 i=1,nziyuu
        do 1100 j=1,nziyuu
          s(i,j)=bc(i)*s(i,j)
          s(j,i)=bc(i)*s(j,i)
1100    continue
        do 1200 i=1,nziyuu
         if( bc(i).lt.1.d-3) then
            s(i,i)=1.d0
            f(i)=0.d0
          end if
1200    continue
c
             p=0.d0
        open(8,file='uvw.out')
        do 2000 l=1,nstep
             p=p+dp
           do 2010 i=1,nziyuu
              df(i)=f(i)*p
           do 2010 j=1,nziyuu
              st(i,j)=s(i,j)
2010       continue
c
c
c 剛性方程式を解く
        call solve(st,d,df,nziyuu,nmax)
c
c 出力
c
       do 3000 j=1,ngamen
        nga=ngame(j)*3-3
         write(*,*) '--------------------------------------------------'
         write(*,*) '節点番号：',ngame(j)
          goto (511,512,513), nhikak
511    write(*,125) ' f=',p,'   PL/EA=',p*ziku/ezz/a
       goto 599
512    continue
        tawa1=p*ziku**3/3.d0/ezz/xi
        tk=10.d0*(1.d0+poiyz)/(12.d0+11.d0*poiyz)
        tawa2=tawa1+p*ziku/gyz/a/tk
       write(*,126) ' f=',p,'    梁理論=',tawa1,'    せん断考慮=',tawa2
       goto 599
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
       do 3010 j=1,nfairu
        nga=ngame(j)*3-3
c        write(8,*) '--------------------------------------------------'
c        write(8,*) '節点番号：',nfair(j)
         write(8,145) p,d(nga+1),d(nga+2),d(nga+3)
3010   continue
c
2000   continue
        close(8)
115      format(a,1pd13.5)
125      format(a,1pd13.5,a,1pd13.5)
126      format(a,1pd13.5,a,1pd13.5,a,1pd13.5)
135      format(a,1p3d13.5)
145      format(1p4d13.5)
           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 代数方程式を解くサブルーチン
c 元は東京大学かどこかの大計のライブラリーか何かを書き換えたもの
      SUBROUTINE SOLVE (ST,DD,DF,N01,nmax)
C
C  GAUSSIAN ELIMINATION METHOD
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ST(nmax,N01),DD(N01),DF(N01)
c
c
C                                 FORWARD ELIMINATION
      DO 100 I=1,N01-1
        INEXT=I+1
        ST(I,I)=1.D0/ST(I,I)
        DF(I)=DF(I)*ST(I,I)
        DO 110 J=INEXT,N01
          ST(I,J)=ST(I,J)*ST(I,I)
          DF(J)=DF(J)-ST(J,I)*DF(I)
          DO 110 K=INEXT,N01
            ST(K,J)=ST(K,J)-ST(K,I)*ST(I,J)
  110   CONTINUE
  100 CONTINUE
C                                 BACKWARD SUBSTITUTION
      DD(N01)=DF(N01)/ST(N01,N01)
      DP=DD(N01)
      DO 120 I=1,N01-1
        K=N01-I
        DD(K)=DF(K)
        KNEXT=K+1
        DO 120 J=KNEXT,N01
          DD(K)=DD(K)-ST(K,J)*DD(J)
  120 CONTINUE
C
      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
           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
        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
