! EUC-koodo
! 8節点24自由度アイソパラメトリック直方体要素を生成するサブルーチンの
! 準備プログラム。
! tyoku*.f 系プログラムを、Melosh要素だけでなくアイソパラメトリック要素
! でも解けるようにすることが目的。
! 05/9/12:ちょっとその前段階の目的として、体積積分とガウス積分と
! その低減積分とMeloshの剛性行列どうしを比較できるようにしたい。
! 06/5/17メモ
!     1)Melosh要素の剛性行列
!     2)アイソパラメトリック要素(原始的な体積積分)
!     3)アイソパラメトリック要素(ガウス積分2積分点)
!     3t)アイソパラメトリック要素(上記の低減積分)
!     3s)アイソパラメトリック要素(上記の選択的低減積分)
!     2s) 2)の選択的低減積分
! で、1),2),3)に有意な差があるようだったら、せん断ロッキングの試算。
! 06/5/18メモ：1)と2)を出力したらまるで合わない。勿論、2)の計算を
! 明らかに間違ってるってことだけど、誰か間違いを捜して…
!
! アイソパラメトリック要素の作り方は、
! 有限要素ハンドブック 7.立体要素、8.アイソパラメトリック要素
! のやり方(表記とか)に従う。
!
! 04/6/1 ようやく準備開始。でも、他にいっぱいやることがあって、
! たぶん、なかなか進まないと思う。どうなることやら。
!
! 私がやろうとしていることを読みとって代わりにやってくれる人が
! 現れたりしてくれると嬉しい……
!
! プログラムを作りかけてる人：後藤文彦
! このプログラムは無保証です。
! このプログラムは改造／再配布して構いません。
! 但し、私的利用に留まらない場合は、もとのプログラムが
! http://www.str.ce.akita-u.ac.jp/~gotou/programoj/
! から取ってきたものと分かるように明示して下さい。
! 改造／再配布後も上記の条件を保持して下さい。
! バグの報告は、
! http://www.str.ce.akita-u.ac.jp/~gotou/png/meeru.png
! にお願いします。
!
        program iso
       implicit real*8 (a-h, o-z)
       dimension xii(8),etai(8),zetai(8),
     & psi38(3,8),xyz83(8,3),yako33(3,3),yakoi33(3,3),
     & di(6,6),d(6,6),yip38(3,8),yiptd83(8,3),yiptdyip88(8,8),
     & si2424(24,24)
! Meloshの剛性行列の計算用
       dimension n8(8),sm2424(24,24),sm(3,3,8,8)
! 
! ******* 材料定数 *********
! x,y,z方向の分割数：nsx,nsy,nsz
       nsx=100; nsy=100; nsz=100
!       nsx=1; nsy=1; nsz=1
! 1要素のx,y,z(横、縦、軸)方向の長さ：a,b,c
       a=1.d0; b=1.d0; c=1.d0
! dV=dxdydz=det[J]dξdηdζのdξ, dη, dζ
       dxi=2.d0/nsx; deta=2.d0/nsy; dzeta=2.d0/nsz
!                  2.d0は積分区間の-1から1までの幅
! dV=dxdydz=det[J]dξdηdζのdξdηdζ
       dxdedz=dxi*deta*dzeta
! 
! ヤング率：exx,eyy,ezz
            exx=1.d0
            eyy=1.d0
            ezz=1.d0
! ポアソン比：poixy,poixz,poiyz
            poixy=1.d0/3.d0
            poixz=1.d0/3.d0
            poiyz=1.d0/3.d0
! せん断弾性係数：gxy,gxz,gyz
            gxy=ezz/2.d0/(1.d0+poixy)
            gyz=gxy
            gzx=gxy
! 有限要素ハンドブックにならってgxzではなくgzxを使っておく。
!
! diに柔性（コンプライアス）行列を入れる
            call lm0(6,6,di)  ! diの初期化
            di(1,1)=1.d0/exx
            di(1,2)=-poixy/exx
            di(1,3)=-poixz/exx
            di(2,1)=di(1,2)
            di(2,2)=1.d0/eyy
            di(2,3)=-poiyz/eyy
            di(3,1)=di(1,3)
            di(3,2)=di(2,3)
            di(3,3)=1.d0/ezz
!
! diの逆行列を取って、応力-ひずみ行列をdに入れる
            call lm0(6,6,d)  !dの初期化
!
            call gyaku(di,d,det)
!
            d(4,4)=gxy
            d(5,5)=gyz
            d(6,6)=gzx
            ! この順番は有限要素ハンドブックに倣っておく。
! 有限要素ハンドブックとMeloshでは、d44,d55,d66の順番が違うので、
! 3つのせん断弾性係数がすべて同じでない場合は、比較の際に注意が必要。
!
! 
! 変位関数の定義に用いた節点番号は、Melosh とは違い、
!      z              5__6
!   ／              ／ ／|
! O →x           1−2   |
! ↓              |  | ／8
! y               3−4
! とする。これに応じて、各節点の無次元化局所座標 xi_i,eta_i,zeta_i を、
! 以下の変数で与えておく。1か-1かの数字で直接 式中に書き込んでもいいと
! 思うけど、式の意味が分かりやすくなるようにするためと、節点番号を
! 変えた場合にも対応できるように変数で与えておく。
!
!
         xii(1)=-1.d0
         xii(2)= 1.d0
         xii(3)=-1.d0
         xii(4)= 1.d0
         xii(5)=-1.d0
         xii(6)= 1.d0
         xii(7)=-1.d0
         xii(8)= 1.d0
         etai(1)= 1.d0
         etai(2)= 1.d0
         etai(3)=-1.d0
         etai(4)=-1.d0
         etai(5)= 1.d0
         etai(6)= 1.d0
         etai(7)=-1.d0
         etai(8)=-1.d0
         zetai(1)=-1.d0
         zetai(2)=-1.d0
         zetai(3)=-1.d0
         zetai(4)=-1.d0
         zetai(5)= 1.d0
         zetai(6)= 1.d0
         zetai(7)= 1.d0
         zetai(8)= 1.d0
!
         call lm0(24,24,si2424) ! si2424の初期化
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
      do zeta=-1,1,dzeta
      do  eta=-1,1,deta
      do   xi=-1,1,dxi
!
!
       do i=1,8
         xyz83(i,1)=a*xii(i)
         xyz83(i,2)=b*etai(i)
         xyz83(i,3)=c*zetai(i)
!
         psi38(1,i)=xii(i)*(1.d0+etai(i)*eta)*(1.d0+zetai(i)*zeta)/8.d0
         psi38(2,i)=etai(i)*(1.d0+xii(i)*xi)*(1.d0+zetai(i)*zeta)/8.d0
         psi38(3,i)=zetai(i)*(1.d0+xii(i)*xi)*(1.d0+etai(i)*eta)/8.d0
       end do
!
! psi38とxyz83をかけてJacobi行列yako33を求める。
         call lmmn(3,8,3,psi38,xyz83,yako33)
! Jacobi行列yako33の逆行列yakoi33と行列式detを求める。
         call gyaku(yako33,yakoi33,det)
!
! yip38 = yakoi33 * psi38 を求める
         call lmmn(3,3,8,yakoi33,psi38,yip38)
!
! k11
         call lm0(8,3,yiptd83)
         call yiptd(yip38,1,1,d(1,1),yiptd83)
         call yiptd(yip38,2,2,d(6,6),yiptd83)
         call yiptd(yip38,3,3,d(5,5),yiptd83)
         call lmmn(8,3,8,yiptd83,yip38,yiptdyip88)
         call si2424kij(1,1,yiptdyip88,det,si2424)
!
! k22
         call lm0(8,3,yiptd83)
         call yiptd(yip38,1,1,d(6,6),yiptd83)
         call yiptd(yip38,2,2,d(2,2),yiptd83)
         call yiptd(yip38,3,3,d(4,4),yiptd83)
         call si2424kij(9,9,yiptdyip88,det,si2424)
!
! k33
         call lm0(8,3,yiptd83)
         call yiptd(yip38,1,1,d(5,5),yiptd83)
         call yiptd(yip38,2,2,d(4,4),yiptd83)
         call yiptd(yip38,3,3,d(3,3),yiptd83)
         call si2424kij(17,17,yiptdyip88,det,si2424)
!
! k12
         call lm0(8,3,yiptd83)
         call yiptd(yip38,1,2,d(1,2),yiptd83)
         call yiptd(yip38,2,1,d(6,6),yiptd83)
         call si2424kij(1,9,yiptdyip88,det,si2424)
!
! k13
         call lm0(8,3,yiptd83)
         call yiptd(yip38,1,3,d(1,3),yiptd83)
         call yiptd(yip38,3,1,d(5,5),yiptd83)
         call si2424kij(1,17,yiptdyip88,det,si2424)
!
! k23
         call lm0(8,3,yiptd83)
         call yiptd(yip38,2,3,d(2,3),yiptd83)
         call yiptd(yip38,3,2,d(4,4),yiptd83)
         call si2424kij(9,17,yiptdyip88,det,si2424)
!
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        end do   ! xi
        end do   ! eta
        end do   ! zeta
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
         do i=1,8
         do j=1,8
          si2424(8+j,i)=si2424(i,8+j)
          si2424(16+j,i)=si2424(i,16+j)
          si2424(16+j,8+i)=si2424(8+i,16+j)
         end do
         end do
!
!
         do i=1,24
         do j=1,24
            si2424(i,j)=si2424(i,j)*dxdedz
         end do
         end do
!
!
         open(7,file='k11i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=1,8),i=1,8)
         close(7)
!
         open(7,file='k12i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=9,16),i=1,8)
         close(7)
!
         open(7,file='k13i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=17,24),i=1,8)
         close(7)
!
         open(7,file='k22i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=9,16),i=9,16)
         close(7)
!
         open(7,file='k22i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=17,24),i=9,16)
         close(7)
!
         open(7,file='k33i.d')
         write(7,'(1p8d10.2)') ((si2424(i,j),j=17,24),i=17,24)
         close(7)
!
! *****************************************************************
! 比較のためにMeloshの剛性行列を計算
          call melesh(a,b,c,d,sm)
! 1要素のsmを、その要素の8節点の節点番号ごとに
! 割り振ってsm2424に書き込む
          call lm0(24,24,sm2424)
        nx=1; ny=1
        i=1; j=1; k=1
         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)
!
        do  l=1,3
        do  m=1,3
        do  ii=1,8
        do  jj=1,8
           i24=(n8(ii)-1)*3
           j24=(n8(jj)-1)*3
            sm2424(i24+l,j24+m)=sm2424(i24+l,j24+m)+sm(l,m,ii,jj)
        end do
        end do
        end do
        end do
!
!
         open(7,file='k11m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=1,8),i=1,8)
         close(7)
!
         open(7,file='k12m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=9,16),i=1,8)
         close(7)
!
         open(7,file='k13m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=17,24),i=1,8)
         close(7)
!
         open(7,file='k22m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=9,16),i=9,16)
         close(7)
!
         open(7,file='k22m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=17,24),i=9,16)
         close(7)
!
         open(7,file='k33m.d')
         write(7,'(1p8d10.2)') ((sm2424(i,j),j=17,24),i=17,24)
         close(7)
!
         end
!
! *****************************************************************
! *****************************************************************
! *****************************************************************
! *****************************************************************
!
!
! l*m行列の初期化
         subroutine lm0(l,m,a)
         implicit real*8(a-h,o-z)
         dimension a(l,m)
         do i=1,l
         do j=1,m
            a(i,j)=0.d0
         end do
         end do
         return
         end
!
!
! 3*3の行列の逆行列と行列式を求める
!
            subroutine gyaku(a,b,det)
             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  i=1,3
             do  j=1,3
              b(i,j)=b(j,i)/det
             end do 
             end do 
             return
             end
!
!
!
! c(l,n)=a(l,m)*b(m,n)のかけ算
         subroutine  lmmn(l,m,n,a,b,c)
         implicit real*8(a-h,o-z)
         dimension a(l,m),b(m,n),c(l,n)
          do i=1,l
          do k=1,n
             c(i,k)=0.d0
          do j=1,m
             c(i,k)=c(i,k)+a(i,j)*b(j,k)
          end do
          end do
          end do
         return
         end
!
!
! 
! yip38^Tのi列目にdを掛け、yiptd83のj列目に代入
         subroutine yiptd(yip38,i,j,d,yiptd83)
         implicit real*8(a-h,o-z)
         dimension yip38(3,8),yiptd83(8,3)
         do k=1,8
         yiptd83(k,j)=yiptd83(k,j)+yip38(i,k)*d
         end do
         return
         end
! 
! 
! si2424に計算したk11,k12,...k33を代入
         subroutine si2424kij(l,m,yiptdyip88,det,si2424)
         implicit real*8(a-h,o-z)
         dimension yiptdyip88(8,8),si2424(24,24)
         do i=l,l+7
         do j=m,m+7
            si2424(i,j)=si2424(i,j)+yiptdyip88(i,j)*det
         end do
         end do
         return
         end
! 
! *****************************************************************
! *****************************************************************
! *****************************************************************
! *****************************************************************
! Meloshの剛性行列の計算(比較用)
        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
