c 集成材梁の上下面に軸方向に溝を掘って鋼リブを挿入して補強した剛性梁の
c 曲げをtyoku2ou.fに半解析させるための入力データを作るプログラム。
c 直方体要素の有限要素プログラム tyoku2ou.f で
c z軸に沿って横たわる片持ち梁の自由端に集中荷重（y方向）を与える問題を
c yz面を対称面とする半解析で解かせるための入力データを作るプログラム。
c このプログラムを実行させると、mage2ou.d というデータができあがるから
c tyoku2ou<mage2ou.d 
c と実行させればよい。
c
c 06/6/12版
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/png/meeru.png
c にお願いします。
c
        program mage2ou
       implicit real*8 (a-h, o-z)
       dimension ijkab(9400),a(3,3),b(3,3),e66a(6,6),e66b(6,6),
     & lkais(100),lryou(100),n8(8),
     & bt1(24,6),bt2(24,6),bt3(24,6),bt4(24,6)
c
c スカイラインの1次元配列のサイズは、同じ自由度数でも要素分割
c に依存するので、自由度数から一意には決められない。
c mage.f の要素分割のやり方だと、軸長（z軸）の分割数がxy断面
c の分割数に比べて多いとスパース部分が増える。
c 仮にスパース部分が皆無だったとすると
c およそ自由度×自由度／２のサイズになるが、こんなサイズは想定外。
c 取り敢えず、12000を最大節点数、その３倍を最大自由度数にしておく。
c 具体的な要素分割例としては、(縦+1)×(横+1)×(軸長+1)が、
c (6+1)*(12+1)*(99+1)=9191,(20+1)*(20+1)*(20+1)=9261 など。
c mage.f で幅（半解析だから幅の半分）×高さ×軸長が
c 1*(8〜20)*(50〜100)ぐらいの範囲のものを解く場合、有効数字３桁
c の収束解が得られる要素分割数は 6分割×10分割×100分割ぐらい。
c 図心が節点上に来るようにするには高さ方向の分割数を偶数にする。
c 幅（の半分）に対して桁高が10倍以上高い場合は、幅の分割数は
c 1 でもいいぐらいかも。
c
c
c 断面横方向の要素分割数：nx 
c 断面縦方向の要素分割数：ny
c 軸方向の要素分割数：nz
          write(*,*) '要素分割数を入力して下さい'
          write(*,*) '（横,縦,軸長の順にコンマ区切りで）'
          write(*,*) '例えば、6,10,100'
          read(*,*) nx,ny,nz
c 鋼材要素の断面横方向の要素分割数：nxa
c 鋼材要素の断面縦方向の要素分割数：nya
          write(*,*) '鋼材要素断面の分割数を入力して下さい'
          write(*,*) '（横,縦の順にコンマ区切りで）'
          write(*,*) '例えば、',1,',',ny/4
          write(*,*) '木材のみにする時は、0,0 '
          write(*,*) '鋼材のみにする時は、',nx,',',ny/2
          write(*,*) '鋼板を上下にぶっとおす時は、',1,',',ny/2
          read(*,*) nxa,nya
c 荷重増分：df
            write(*,*) '荷重増分を入力して下さい'
            read(*,*) df
c             df=1.d0
c 荷重増分ステップ数：nstep
c            write(*,*) '荷重増分ステップ数を入力して下さい'
c            read(*,*) nstep
                       nstep=1
c 梁断面の幅：yoko, 桁高：tate
            yoko=30.d-3
            tate=120.d-3
c 鋼板の断面の横の長さ:yokoa0, 縦の長さ:tatea0
          if( (nxa.eq.nx).and.(nya.eq.ny/2) ) then
c         もし全部鋼材だったら
            yokoa0=yoko
            tatea0=tate/2
          else
c         鋼材が全部ではないとき
            yokoa0=4.5d-3
            tatea0=40.d-3
          end if
c 鋼材要素の断面の横の長さ:yokoa, 縦の長さ:tatea
          if( (nxa.eq.0).and.(nya.eq.0) ) then
c         もし木材だけだったら
            yokoa0=0.d0
            tatea0=0.d0
            yokoa=0.d0
            tatea=0.d0
          else
c         鋼材があるときは
            yokoa=yokoa0/float(nxa)
            tatea=tatea0/float(nya)
          end if
c 木材要素(分類番号4）の断面の横の長さ:yokob, 縦の長さ:tateb
          if( (nxa.eq.nx).and.(nya.eq.ny/2) ) then
c         もし全部鋼材だったら
            yokob=0.d0
            tateb=0.d0
          else
c         鋼材が全部ではないとき
            if(nya.eq.ny/2) then
c           もし鋼板が上から下までぶっとおっていたら
             yokob=(yoko-yokoa0)/(float(nx)-float(nxa))
             tateb=0.d0
            else
             yokob=(yoko-yokoa0)/(float(nx)-float(nxa))
             tateb=(tate-2.d0*tatea0)/(float(ny)-2*float(nya))
            end if
          end if
c 鋼板が底面に貼ってある（左右にぶっとおっている）ようなものは、
c aを木材,bを鋼材に置き換えて90゜回転させて解く。
c
c 梁の軸長：ziku
            ziku=1080.d-3
            ziku1=ziku/float(nz)
            write(*,*) 'yokoa=',yokoa
            write(*,*) 'tatea=',tatea
            write(*,*) 'yokob=',yokob
            write(*,*) 'tateb=',tateb
            write(*,*) 'ziku1=',ziku1
c
c 要素番号(i,j,k)の要素の分類番号
c yokoa*tateaの鋼：1
c yokob*tateaの木：2
c yokoa*tatebの木：3
c yokob*tatebの木：4
          call bunrui(nx,ny,nz,nxa,nya,ijkab)
c
c
c ********** 境界条件 *************
c 
c 拘束箇所の数：kkasyo
c 同じ節点を定義しなおした場合は、後から定義した方が有効
c 一箇所目で定義した場所を二箇所目で定義しなおすこともできる
c 固定端の面が1箇所目
c 固定端の中立軸が2箇所目
c 固定端の対称軸が3〜3+ny箇所目
c 対称面が3+ny+1〜3+ny+(ny+1)*nz箇所目
c 中立軸対称点が 3+ny+(ny+1)*nz+1箇所目
            kkasyo=3+ny+(ny+1)*nz+1
c 1箇所目の拘束箇所における拘束開始節点番号:kkais1
            kkais1=1
c 1箇所目の拘束箇所における拘束終了節点番号:kryou1
            kryou1=(nx+1)*(ny+1)
c tyoku.fでは拘束開始点から拘束終了点まで節点番号順に全て
c 同じ拘束条件を与える。要素番号の順番は、最前面最上段左端の節点を
c 1として、右(x)方向に進み、右端まできたらy方向へ1段下の段の左端から
c 右へすすみ、真下の右端までいったら、z方向へ1面奥の面の左上から、
c また同じように1段づつ右へ進み、最下段の右端までいったら、もう一つ
c 奥の面へと進んでいき、一番奥の面の最下段の右端が最終要素番号
c 1箇所目の拘束箇所における拘束条件
c u,v,w 方向の変位をそれぞれ拘束するか開放するか(0/1)
c まず、端部のz方向変位のみをすべて拘束する
             bu1=1.d0
             bv1=1.d0
             bw1=0.d0
c でも、xy平面を移動できてしまうので、中立軸のy方向変位を拘束する
             kkais2=(nx+1)*int(float(ny)/2.d0)+1
             kryou2=kkais2+nx
             bu2=1.d0
             bv2=0.d0
             bw2=0.d0
c あと、対称面(yz面)のx方向変位を拘束する必要があるが、これは、
c 量が多いので下のデータ出力時にdo ループで直接 書く。
c
c
c *********** 載荷条件 ************
c
c 載荷箇所の数：lkasyo
c 載荷箇所は、同じ節点の荷重を後から定義しなおしても構わない。
c 今回は、2要素の重なる節点を1箇所目、
c 四隅の1要素だけの節点を2,3,4,5箇所目
c 4要素の重なる節点を6箇所目として与える
c（縦横の分割数が共に2以上のときのみ）。
c 
           if((nx.gt.1).and.(ny.gt.1)) then
            lkasyo=5+ny-1
           else
            lkasyo=5
           endif
c 自由端に一様な引張を与えるのだけど、節点をn個の要素で
c 共有している節点には n倍の節点力を与える。
c 大きさ 1 の荷重をのべ何個の節点に振り分けるか：nf
             nf=4*nx*ny
             fn=float(nf)
c まず最初に、自由端の全ての節点に fy=2.d0 を与えてしまう。
c 1箇所目の載荷箇所における載荷開始節点番号:lkais(1)
c 自由端の節点番号の一番若い節点（左上のかど）
            lkais(1)=(nx+1)*(ny+1)*nz+1
c 1箇所目の載荷箇所における載荷終了節点番号:lryou(1)
            lryou(1)=(nx+1)*(ny+1)*(nz+1)
c tyoku.fでは載荷開始点から載荷終了点まで節点番号順に全て
c 同じ載荷条件を与える。要素番号の順番は、最前面最上段左端の節点を
c 1として、右(x)方向に進み、右端まできたらy方向へ1段下の段の左端から
c 右へすすみ、真下の右端までいったら、z方向へ1面奥の面の左上から、
c また同じように1段づつ右へ進み、最下段の右端までいったら、もう一つ
c 奥の面へと進んでいき、一番奥の面の最下段の右端が最終要素番号
c 1箇所目の載荷箇所における載荷条件
c fx,fy,fz 方向のそれぞれに与える荷重
             fx2=0.d0
             fy2=2.d0/fn
             fz2=0.d0
c
             lkais(2)=lkais(1)
             lryou(2)=lkais(1)
             lkais(3)=lkais(2)+nx
             lryou(3)=lkais(3)
             lkais(4)=lkais(1)+(nx+1)*ny
             lryou(4)=lkais(4)
             lkais(5)=lryou(1)
             lryou(5)=lkais(5)
              fx1=0.d0
              fy1=1.d0/fn
              fz1=0.d0
c
c 次に、自由端の周を除くすべての節点に fz=4.d0 を与えなおす。
           if((nx.gt.1).and.(ny.gt.1)) then
            do 50 i=1,ny-1
               lkais(5+i)=lkais(1)+(nx+1)*i+1
               lryou(5+i)=lkais(1)+(nx+1)*i+nx-1
50          continue
             fx4=0.d0
             fy4=4.d0/fn
             fz4=0.d0
            endif
c
c ******* 材料定数 *********
c 鋼材の材料定数
c ヤング率：exxa,eyya,ezza
            exxa=200.d9
            eyya=200.d9
            ezza=200.d9
c ポアソン比：poixya,poixza,poiyza
            poixya=1.d0/3.d0
            poixza=1.d0/3.d0
            poiyza=1.d0/3.d0
c せん断弾性係数：gxya,gxza,gyz
            gxya=ezza/2.d0/(1.d0+poixya)
            gxza=gxya
            gyza=gxya
c
c aに柔性（コンプライアス）行列を入れる
            a(1,1)=1.d0/exxa
            a(1,2)=-poixya/exxa
            a(1,3)=-poixza/exxa
            a(2,1)=a(1,2)
            a(2,2)=1.d0/eyya
            a(2,3)=-poiyza/eyya
            a(3,1)=a(1,3)
            a(3,2)=a(2,3)
            a(3,3)=1.d0/ezza
c
c 木材の材料定数
c ヤング率：exxb,eyyb,ezzb
            exxb=0.5d9
            eyyb=1.d9
            ezzb=10.d9
c ポアソン比：poixyb,poixzb,poiyzb
            poixyb=0.2d0
            poixzb=0.025d0
            poiyzb=0.03d0
c せん断弾性係数：gxyb,gxzb,gyzb
            gxyb=0.7d9
            gxzb=0.7d9
            gyzb=0.7d9
c
c bに柔性（コンプライアス）行列を入れる
            b(1,1)=1.d0/exxb
            b(1,2)=-poixyb/exxb
            b(1,3)=-poixzb/exxb
            b(2,1)=b(1,2)
            b(2,2)=1.d0/eyyb
            b(2,3)=-poiyzb/eyyb
            b(3,1)=b(1,3)
            b(3,2)=b(2,3)
            b(3,3)=1.d0/ezzb
c
c 鋼材と木材の柔性行列の逆行列と、せん断弾性係数をそれぞれ
c 66a,e66bに代入し、それぞれの応力−ひずみ行列を作る。
            call gyaku(a,gxya,gxza,gyza,e66a)
            call gyaku(b,gxyb,gxzb,gyzb,e66b)
c
c 要素ごとのBマトリックスを求める
c nは以下の分類番号
c yokoa*tateaの鋼：n=1
c yokob*tateaの木：n=2
c yokoa*tatebの木：n=3
c yokob*tatebの木：n=4
         call ondo(yokoa,tatea,ziku1,e66a,1,bt1)
         call ondo(yokob,tatea,ziku1,e66b,2,bt2)
         call ondo(yokoa,tateb,ziku1,e66b,3,bt3)
         call ondo(yokob,tateb,ziku1,e66b,4,bt4)
c
c
c tyokuon.fで出力された節点変位（fort.12に出力）に、ouryoku.fで
c 歪-変位行列と応力-歪行列をかけて各要素の応力を求める時のために、
c 要素分割数、要素ごとの分類番号、分類番号ごとの歪-変位行列を
c フォーマットなしでfort.11に書き出しておく。
       write(11) nx,ny,nz,ijkab,e66a,e66b,bt1,bt2,bt3,bt4
       rewind(11)
c
c ***** 出力の設定 *****
c 比較のために出力したい梁−柱理論値：nhikak
c nhikak=1: w 方向に引っ張った時の先端の伸びを出力
c nhikak=2: y 方向に曲げた時の先端のたわみを出力
c nhikak=3: x 方向に曲げた時の先端のたわみを出力
                  nhikak=2
c
c 変位を画面出力したい節点の数：ngamen
c 例えば、先端の節点すべてだったら：
            ngamen=(nx+1)*(ny+1)
c 1つ目の出力したい節点番号：ngame1
            ngame1=(nx+1)*(ny+1)*nz+1
c
c 変位をファイル出力したい節点の数：nfairu
            nfairu=1
c 1つ目の出力したい節点番号：nfair1
            nfair1=(nx+1)*(ny+1)*nz
     &            +(nx+1)*int(float(ny)/2.d0)
     &            +1
c
c
c
c tyoku.f 用の入力データを kata.d というファイルに出力する
c
         open(7,file='mage2ou.d')
         write(7,1000) nx
         write(7,1000) ny
         write(7,1000) nz
          nyou=nx*ny*nz
         write(7,1000) nxa
         write(7,1000) nya
          do 70 i=1,nyou
         write(7,1000) ijkab(i)
70        continue
         write(7,500) yokoa
         write(7,500) tatea
         write(7,500) yokob
         write(7,500) tateb
         write(7,500) ziku
         write(7,1000) kkasyo
         write(7,1000) kkais1
         write(7,1000) kryou1
         write(7,500) bu1
         write(7,500) bv1
         write(7,500) bw1
         write(7,1000) kkais2
         write(7,1000) kryou2
         write(7,500) bu2
         write(7,500) bv2
         write(7,500) bw2
c 固定端の対称面(y軸)のx,z方向変位を固定
        do 30 j=0,ny
         write(7,1000) (nx+1)*j+1
         write(7,1000) (nx+1)*j+1
         write(7,500) 0.d0
         write(7,500) 1.d0
         write(7,500) 0.d0
30      continue
c 対称面(yz面)のx方向変位を固定
        do 40 k=1,nz
        do 40 j=0,ny
         write(7,1000) (nx+1)*(ny+1)*k+(nx+1)*j+1
         write(7,1000) (nx+1)*(ny+1)*k+(nx+1)*j+1
         write(7,500) 0.d0
         write(7,500) 1.d0
         write(7,500) 1.d0
40      continue
c 固定端中立軸の対称点(x=0の点）を固定
         write(7,1000) kkais2
         write(7,1000) kkais2
         write(7,500) 0.d0
         write(7,500) 0.d0
         write(7,500) 0.d0
c
         write(7,1000) lkasyo
         write(7,1000) lkais(1)
         write(7,1000) lryou(1)
         write(7,500) fx2
         write(7,500) fy2
         write(7,500) fz2
         write(7,1000) lkais(2)
         write(7,1000) lryou(2)
         write(7,500) fx1
         write(7,500) fy1
         write(7,500) fz1
         write(7,1000) lkais(3)
         write(7,1000) lryou(3)
         write(7,500) fx1
         write(7,500) fy1
         write(7,500) fz1
         write(7,1000) lkais(4)
         write(7,1000) lryou(4)
         write(7,500) fx1
         write(7,500) fy1
         write(7,500) fz1
         write(7,1000) lkais(5)
         write(7,1000) lryou(5)
         write(7,500) fx1
         write(7,500) fy1
         write(7,500) fz1
           if((nx.gt.1).and.(ny.gt.1)) then
            do 20 i=1,ny-1
         write(7,1000) lkais(5+i)
         write(7,1000) lryou(5+i)
         write(7,500) fx4
         write(7,500) fy4
         write(7,500) fz4
20          continue
           endif
c
        do 10 i=1,3
        do 10 j=1,3
         write(7,900) a(i,j)
10      continue
         write(7,500) gxya
         write(7,500) gxza
         write(7,500) gyza
c
        do 11 i=1,3
        do 11 j=1,3
         write(7,900) b(i,j)
11      continue
         write(7,500) gxyb
         write(7,500) gxzb
         write(7,500) gyzb
c
         write(7,500) df
         write(7,1000) nhikak
c 端部の図心変位のみを出す場合
         write(7,1000) nfairu
         write(7,1000) nfair1
c 端部の全節点の変位を出す場合
c         write(7,1000) ngamen
c          do 21 i=lkais(1),lryou(1)
c           write(7,1000) i
c21        continue
         write(7,1000) nfairu
         write(7,1000) nfair1
         write(7,1000) nstep
          close(7)
c 小数点以下5桁の実数
500      format(1pd13.5)
c 小数点以下9桁の実数
900      format(1pd17.9)
c 10桁までの整数
1000     format(i10)
         stop
         end
c
c
c 要素番号(i,j,k)の要素の分類番号を決めるサブルーチン
c yokoa*tateaの鋼：1
c yokob*tateaの木：2
c yokoa*tatebの木：3
c yokob*tatebの木：4
c
         subroutine bunrui(nx,ny,nz,nxa,nya,ijkab)
         implicit real*8 (a-h, o-z)
         dimension ijkab(9400)
c もしnxa=0,nya=0（鋼材なし）だったら100に飛ぶ
         if( (nxa.eq.0).and.(nya.eq.0)) goto 100
c
c
c 左上の鋼板要素：横nxa*縦nya分割
            do i=1,nxa
            do j=1,nya
            do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=1
            end do
            end do
            end do
c 左下の鋼板要素：横nxa*縦nya分割
            do i=1,nxa
               ja=ny-nya+1
            do j=ja, ny
            do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=1
            end do
            end do
            end do
c          もしnxa=nx,nya=ny/2(全部鋼材)だったら、ここで戻る
          if( (nxa.eq.nx).and.(nya.eq.ny/2) ) return
c
c 上鋼板右隣の木材要素：横(nx-nxa)*縦(nya)分割
            do i=nxa+1,nx
            do j=1,nya
            do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=2
            end do
            end do
            end do
c 下鋼板右隣の木材要素：横(nx-nxa)*縦(nya)分割
           do i=nxa+1,nx
              ja=ny-nya+1
           do j=ja, ny
           do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=2
            end do
            end do
            end do
c もし鋼板が上から下までぶっとおっていたら、ここで戻る
         if(nya.eq.ny/2) return
c
c 鋼板に挟まれた左端の木材要素：横nx*縦(ny-2*nya)分割
           do i=1,nxa
           do j=nya+1,ny-nya
           do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=3
           end do
           end do
           end do
c
100        continue
c 残りの木材要素：横(nx-nxa)*縦(ny-nya)分割
           do i=nxa+1,nx
           do j=nya+1,ny-nya
           do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=4
           end do
           end do
           end do
          return
          end
c
c
c **********************************************************
c ひずみ−応力行列aの逆数を取って、応力−ひずみ行列を求め、
c それとせん断弾性係数gxy,gxz,gyzとをeに代入して、
c 応力−ひずみ行列を作る。
c
            subroutine gyaku(a,gxy,gxz,gyz,e)
             implicit real*8(a-h, o-z)
             dimension a(3,3),b(3,3),e(6,6)
              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
              e(i,j)=b(j,i)/det
              e(i,3+j)=0.d0
              e(3+i,j)=0.d0
              e(3+i,3+j)=0.d0
             end do
             end do
              e(4,4)=gxy
              e(5,5)=gxz
              e(6,6)=gyz
             return
             end
c
c
c 
c 歪−変位行列（いわゆるBマトリクス）の1要素の体積積分を求め、
c それの転置を bt とする
c 変位関数の定義に用いた変位d1-d8の番号は、Melosh とは違い、
c      z              5__6
c   ／              ／ ／|
c O →x           1−2   |
c ↓              |  | ／8
c y               3−4
c だから、荷重ベクトルの重ね合わせで使っているn8の添え字が、
c tyoku2.f とは違っている。
         subroutine ondo(yoko1,tate1,ziku1,e66,nab,bt)
         implicit real*8 (a-h, o-z)
         dimension bt(24,6),e66(6,6)
          do i=1,24
          do j=1,6
             bt(i,j)=0.d0
          end do
          end do
c
           a=yoko1/2.d0
           b=tate1/2.d0
           c=ziku1/2.d0
          bt( 1,1)=-b*c
          bt( 4,1)= b*c
          bt( 7,1)=-b*c
          bt(10,1)= b*c
          bt(13,1)=-b*c
          bt(16,1)= b*c
          bt(19,1)=-b*c
          bt(22,1)= b*c
c
          bt( 2,2)=-a*c
          bt( 5,2)=-a*c
          bt( 8,2)=+a*c
          bt(11,2)=+a*c
          bt(14,2)=-a*c
          bt(17,2)=-a*c
          bt(20,2)=+a*c
          bt(23,2)=+a*c
c
          bt( 3,3)=-a*b
          bt( 6,3)=-a*b
          bt( 9,3)=-a*b
          bt(12,3)=-a*b
          bt(15,3)= a*b
          bt(18,3)= a*b
          bt(21,3)= a*b
          bt(24,3)= a*b
c
c 温度によるせん断歪がないなら、tyokuon.fで節点変位を求める上では、
c bt(i,j)のj>4の成分も必要は必要ないだろうが、
c 鋼と木の熱膨張が異なることによってせん断歪が生じるなら、
c ouryoku.fで応力を求める時のために以下の成分が必要になる。
c
          bt( 1,4)= a*c
          bt( 2,4)=-b*c
          bt( 4,4)= a*c
          bt( 5,4)= b*c
          bt( 7,4)=-a*c
          bt( 8,4)=-b*c
          bt(10,4)=-a*c
          bt(11,4)= b*c
          bt(13,4)= a*c
          bt(14,4)=-b*c
          bt(16,4)= a*c
          bt(17,4)= b*c
          bt(19,4)=-a*c
          bt(20,4)=-b*c
          bt(22,4)=-a*c
          bt(23,4)= b*c
c
          bt( 1,5)=-a*b
          bt( 3,5)=-b*c
          bt( 4,5)=-a*b
          bt( 6,5)= b*c
          bt( 7,5)=-a*b
          bt( 9,5)=-b*c
          bt(10,5)=-a*b
          bt(12,5)= b*c
          bt(13,5)= a*b
          bt(15,5)=-b*c
          bt(16,5)= a*b
          bt(18,5)= b*c
          bt(19,5)= a*b
          bt(21,5)=-b*c
          bt(22,5)= a*b
          bt(24,5)= b*c
c
          bt( 2,6)=-a*b
          bt( 3,6)= a*c
          bt( 5,6)=-a*b
          bt( 6,6)= a*c
          bt( 8,6)=-a*b
          bt( 9,6)=-a*c
          bt(11,6)=-a*b
          bt(12,6)=-a*c
          bt(14,6)= a*b
          bt(15,6)= a*c
          bt(17,6)= a*b
          bt(18,6)= a*c
          bt(20,6)= a*b
          bt(21,6)=-a*c
          bt(23,6)= a*b
          bt(24,6)=-a*c
c
         do i=1,24
         do j=1,6
          bt(i,j)=bt(i,j)/a/b/c/8.d0
         end do
         end do
c
       return
       end
c
c
c
