c hippon.f を改造中．
c
c
c 集成材梁の上下面に軸方向に溝を掘って鋼リブを挿入して補強した剛性梁に
c 温度（と含水による膨張）による荷重を与えて（右下だけを）tyokuon.f に
c 1/4解析させるための入力データを作るプログラム。
c 計算例は
c 後藤文彦，高橋明洋，薄木征三，佐々木貴信： 
c 熱膨張・乾燥収縮を受ける鋼板挿入集成 材梁の有限要素解析， 
c 第2回木橋技術に関するシンポジウム論文報告集， pp. 97-100, (2003).
c とかだけど、鋼板・集成材境界面の要素の変形が大きすぎるので、
c あまりいい計算例ではない。
c
c   O___________x(上下対称面）
c    LLLLLL
c    LLLLLL
c    LLLLLL
c    #}LLLL         Lが集成材要素（場所によってサイズは3種類）
c    #}LLLL         #が鋼材要素（サイズは1種類）
c   |               }が接着剤要素（サイズは1種類）
c   y（左右対称面）
c
c 直方体要素の有限要素プログラム tyokuon.f で
c z軸に沿って横たわる片持ち梁の自由端に集中荷重（y方向）を与える問題を
c yz面を対称面とする半解析で解かせるための入力データを作るプログラム。
c このプログラムを実行させると、hippon.d というデータができあがるから
c tyokuon<hippon.d (hogeにファイル出力するなら tyokuon<hippon.d>hoge ) 
c と実行させればよい。
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 03/10/9版
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 hippou
       implicit real*8 (a-h, o-z)
       dimension ijkab(9400),a(3,3),b(3,3),c(3,3),
     & e66a(6,6),e66b(6,6),e66c(6,6),
     & f(9400*3),fab(6,24),n8(8),
     & bt1(24,6),bt2(24,6),bt3(24,6),bt4(24,6),bt5(24,6),bt6(24,6)
c
c スカイラインの1次元配列のサイズは、同じ自由度数でも要素分割
c に依存するので、自由度数から一意には決められない。
c mage.f の要素分割のやり方だと、軸長（z軸）の分割数がxy断面
c の分割数に比べて多いとスパース部分が増える。
c 仮にスパース部分が皆無だったとすると
c およそ自由度×自由度／２のサイズになるが、こんなサイズは想定外。
c 取り敢えず、9400を最大節点数、その３倍を最大自由度数にしておく。
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 断面横方向の要素分割数：nx 
c 断面縦方向の要素分割数：ny
c 軸方向の要素分割数：nz
          write(*,*) '要素分割数を入力して下さい'
          write(*,*) '（横,縦,軸長の順にコンマ区切りで）'
          write(*,*) '例えば、20,4,2'
          read(*,*) nx,ny,nz
c 接着剤要素の断面横方向の要素分割数：nxc
c 接着剤要素の断面縦方向の要素分割数：nyc
          write(*,*) ' 接着剤要素断面の分割数を入力して下さい'
          write(*,*) ' （横,縦の順にコンマ区切りで）'
c          write(*,*) ' 例えば、',2,',',ny/2
          write(*,*) '鋼板を上下にぶっとおす時の縦の分割数は',ny
          read(*,*) nxc,nyc
c 鋼材要素の断面横方向の要素分割数：nxa
c 鋼材要素の断面縦方向の要素分割数：nya
           nya=nyc
          write(*,*) '鋼材要素の幅方向の分割数を入力して下さい'
          write(*,*) ' 標準は 1'
          read(*,*) nxa
c 荷重増分：df
c            write(*,*) '荷重増分を入力して下さい'
c            read(*,*) df
             df=1.d0
c 荷重増分ステップ数：nstep
c            write(*,*) '荷重増分ステップ数を入力して下さい'
c            read(*,*) nstep
                       nstep=1
c 温度変化：dt
            write(*,*) '温度変化を（℃で）入力して下さい'
            read(*,*) dt
c 木材含水率変化：dh
           write(*,*) '含水率変化を（％で）入力して下さい'
           read(*,*) dh
c 梁断面の幅：yoko, 桁高：tate
            yoko=100.d-3
            tate=400.d-3
c 鋼板の断面の横の長さ:yokoa0, 縦の長さ:tatea0
c 接着剤層の横の長さ:yokoc0, 縦の長さ:tatec0
          if( (nxa.eq.nx).and.(nya.eq.ny) ) then
c         もし全部鋼材だったら
            yokoa0=yoko
            tatea0=tate
          else
c         鋼材が全部ではないとき
            yokoa0=5.d-3
            tatea0=400.d-3
            yokoc0=1.d-3
            tatec0=tatea0
          end if
c 鋼材要素の断面の横の長さ:yokoa, 縦の長さ:tatea
c 接着剤要素の断面の横の長さ:yokoc, 縦の長さ:tatec
          if( (nxa.eq.0).and.(nya.eq.0) ) then
c         もし木材だけだったら
            yokoa0=0.d0
            tatea0=0.d0
            yokoa=0.d0
            tatea=0.d0
            yokoc0=0.d0
            tatec0=0.d0
            yokoc=0.d0
            tatec=0.d0
          else
c         鋼材があるときは
            yokoa=yokoa0/nxa
            tatea=tatea0/nya
            yokoc=yokoc0/nxc
            tatec=tatec0/nyc
          end if
c 木材要素(分類番号4）の断面の横の長さ:yokob, 縦の長さ:tateb
          if( (nxa.eq.nx).and.(nya.eq.ny) ) then
c         もし全部鋼材だったら
            yokob=0.d0
            tateb=0.d0
            yokoc=0.d0
            tatec=0.d0
          else
c         鋼材が全部ではないとき
            if(nya.eq.ny) then
c           もし鋼板が上から下までぶっとおっていたら
             yokob=(yoko-yokoa0-yokoc0)/(nx-nxa-nxc)
             tateb=0.d0
            else
             yokob=(yoko-yokoa0-yokoc0)/(nx-nxa-nxc)
             tateb=(tate-tatea0)/(ny-nya)
            end if
          end if
c 鋼板が底面に貼ってある（左右にぶっとおっている）ようなものは、
c aを木材,bを鋼材に置き換えて90゜回転させて解く。
c
c 梁の軸長：ziku
            ziku=10.d0
            ziku1=ziku/float(nz)
            write(*,*) 'yokoa=',yokoa
            write(*,*) 'tatea=',tatea
            write(*,*) 'yokob=',yokob
            write(*,*) 'tateb=',tateb
            write(*,*) 'yokoc=',yokoc
            write(*,*) 'tatec=',tatec
            write(*,*) 'ziku1=',ziku1
          call bunrui(nx,ny,nz,nxa,nya,nxc,ijkab)
c
c ********** 境界条件 *************
c
c 拘束箇所の数：nkasyo
c 同じ節点を定義しなおした場合は、後から定義した方が有効
c 一箇所目で定義した場所を二箇所目で定義しなおすこともできる
c 固定端の面（z=0面）が1箇所
c 固定端面の上辺（x軸）が2箇所目
c 左右対称面（y=0面）が3～3+nz箇所目
c 固定端面の左端（y軸）が3+nz～3+nz+ny箇所目
c 上下対称面（x=0面）が3+nz+ny+1～3+nz+ny+(ny+1)*nz箇所目
c 対称軸（z軸）が3+nz+ny+(ny+1)*nz+1～3+nz+ny+(ny+1)*nz+nz箇所目
c 対称点（原点）が 3+nz+ny+(ny+1)*nz+nz+1箇所目
            kkasyo=3+nz+ny+(ny+1)*nz+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 固定端の上下対称線（上辺）のyz方向を拘束する
             kkais2=1
             kryou2=nx+1
             bu2=1.d0
             bv2=0.d0
             bw2=0.d0
c 上下対称面、左右対称面の拘束箇所はいっぱいあるので、
c 下のデータ出力のところで、doループで直接 書き出す
c
c 因みに、固定端の全ての節点でxy方向の変位を拘束してしまうと、
c （つまり、固定端で材料がちぢめないような境界条件だと）
c w=LP/EAの解からは、ちょっとずれる（要素数を増やしても）。
c
c
c
c ******* 材料定数 *********
c 鋼材の材料定数
c ヤング率：exxa,eyya,ezza
            exxa=200.d9
            eyya=exxa
            ezza=exxa
c ポアソン比：poixya,poixza,poiyza
            poixya=1.d0/3.d0
            poixza=poixya
            poiyza=poixya
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 鋼材の線膨張率：aa
            aa=12.d-6
c 鋼材の温度歪：epa
            epa=aa*dt
c
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.03d0
            gxzb=0.6d0
            gyzb=0.7d0
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 木材の軸（繊維）方向の線膨張率：tbz
            tbz=3.5d-6
c 木材の軸直角方向の線膨張率：tbx,tby
            tbx=30.d-6
            tby=30.d-6
c 参考：産業調査会『木材活用事典』（事典出版センター）
c
c 木材の軸（繊維）方向の収縮率（含水率1%当たり）：hbz
            hbz=0.015d-2
c 木材の軸直角方向の収縮率（含水率1%当たり）：hbx,hby
            hbx=0.2d-2
            hby=0.2d-2
c 参考：木材活用辞典,(株)産業調査会辞典出版センター(1994)
c http://www.cfqlcs.go.jp/technical_information/tree_and_life/wl38.htm
c
c 温度と湿度の変化による木材の歪：epbx,epby,epbz
            epbx=tbx*dt + hbx*dh
            epby=tby*dt + hby*dh
            epbz=tbz*dt + hbz*dh
c
c           write(*,*) 'epa,epbx,epby,epbz=',epa,epbx,epby,epbz
c
c 接着剤の材料定数
c 木質複合建築構造技術の開発平成１３年度報告書－構造分科会－
c 国土交通省国土技術政策総合研究所、独立行政法人建築研究所、
c （財）日本建築センター, p.101-p.108
c ２．６　接着耐久性に関する検討
c ２．６．３　実験および解析
c （１）ペンギンセメント#1031の木材／ブラスと処理鋼板接着性
c および樹脂硬化物の物性について
c
c ヤング率：exxc,eyyc,ezzc
            exxc=2.4d9
            eyyc=exxc
            ezzc=exxc
c ポアソン比：poixya,poixza,poiyza
            poixyc=0.41
            poixzc=poixyc
            poiyzc=poixyc
c せん断弾性係数：gxya,gxza,gyz
            gxyc=ezzc/2.d0/(1.d0+poixyc)
            gxzc=gxyc
            gyzc=gxyc
c
c aに柔性（コンプライアス）行列を入れる
            c(1,1)=1.d0/exxc
            c(1,2)=-poixyc/exxc
            c(1,3)=-poixzc/exxc
            c(2,1)=c(1,2)
            c(2,2)=1.d0/eyyc
            c(2,3)=-poiyzc/eyyc
            c(3,1)=c(1,3)
            c(3,2)=c(2,3)
            c(3,3)=1.d0/ezzc
c
c サンスター技研株式会社
c ケミカル事業部新規開発グループマーケティング部
c 【＃1090 の特徴】
c http://www.sunstar-engineering.com/sunstar-engineering/
c chemical/pdf/secondary.pdf
c
c 接着剤の線膨張率：ac
            ac=100.d-6
c 接着剤の温度歪：epc
            epc=ac*dt
c
c 鋼材、接着剤、木材の柔性行列の逆行列と、せん断弾性係数をそれぞれ
c 66a,e66c,e66bに代入し、それぞれの応力－ひずみ行列を作る。
            call gyaku(a,gxya,gxza,gyza,e66a)
            call gyaku(b,gxyb,gxzb,gyzb,e66b)
            call gyaku(c,gxyc,gxzc,gyzc,e66c)
c
c 要素ごとの節点温度（湿度）荷重ベクトル：fab(n,24)を求める
c nは以下の分類番号
c yokoa*tateaの鋼材：  1
c yokoc*tateaの接着剤：2
c yokob*tateaの集成材：3
c yokoa*tatebの集成材：4
c yokoc*tatebの集成材：5
c yokob*tatebの集成材：6
         call ondo(yokoa,tatea,ziku1,e66a, epa, epa, epa,fab,1,bt1)
         call ondo(yokoc,tatea,ziku1,e66c, epc, epc, epc,fab,2,bt2)
         call ondo(yokob,tatea,ziku1,e66b,epbx,epby,epbz,fab,3,bt3)
         call ondo(yokoa,tateb,ziku1,e66b,epbx,epby,epbz,fab,4,bt4)
         call ondo(yokoc,tateb,ziku1,e66b,epbx,epby,epbz,fab,5,bt4)
         call ondo(yokob,tateb,ziku1,e66b,epbx,epby,epbz,fab,6,bt6)
c
       open(9,file='ondo.d')
       write(9,*) '分類番号1～4の要素の節点温度荷重'
       write(9,*) 'fx'
       do i=1,22,3
      write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
       end do
       write(9,*) 'fy'
       do i=2,23,3
      write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
       end do
       write(9,*) 'fz'
       do i=3,24,3
      write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
       end do
       close(9)
9      format(1p6d11.3)
c
c tyokuon.fで出力された節点変位（fort.12に出力）に、ouryogu.fで
c 歪-変位行列と応力-歪行列をかけて各要素の応力を求める時のために、
c 要素分割数、要素ごとの分類番号、分類番号ごとの歪-変位行列を
c フォーマットなしでfort.11に書き出しておく。
      write(11) nx,ny,nz,ijkab,e66a,e66b,e66c,bt1,bt2,bt3,bt4,bt5,bt6
      rewind(11)
c
c 要素ごとの節点温度荷重を（要素の8節点の節点番号ごとに割り振って）
c 重ね合わせて全体の温度荷重ベクトル：f
c を求める
c 全節点数：nset
          nset=(nx+1)*(ny+1)*(nz+1)
c 全自由度数：nziyuu
          nziyuu=nset*3
         do 71 i=1,nziyuu
            f(i)=0.d0
71       continue
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 1001 m=1,3
        do 1001 jj=1,8
           j24=(n8(jj)-1)*3
c            s(i24+l,j24+m)=s(i24+l,j24+m)+sm(l,m,ii,jj)
             jj1=jj-1
            f(j24+m)=f(j24+m)+fab(nab,3*jj1+m)
1001    continue
c
c
c ***** 出力の設定 *****
c 比較のために出力したい梁－柱理論値：nhikak
c nhikak=1: w 方向に引っ張った時の先端の伸びを出力
c nhikak=2: y 方向に曲げた時の先端のたわみを出力
c nhikak=3: x 方向に曲げた時の先端のたわみを出力
                  nhikak=2
c
c
c 変位を画面出力したい節点の数：ngamen
            ngamen=1
c 1つ目の出力したい節点番号：ngame1
            ngame1=(nx+1)*(ny+1)*nz
     &            +1
c
c 変位をファイル出力したい節点の数：nfairu
c 例えば、先端の節点すべてだったら：
            nfairu=(nx+1)*(ny+1)
c 1つ目の出力したい節点番号：nfair1
            nfair1=(nx+1)*(ny+1)*nz+1
c
c
c hippon.f 用の入力データを hippon.d というファイルに出力する
c
         open(7,file='hippon3.d')
         write(7,1000) nx
         write(7,1000) ny
         write(7,1000) nz
c 全要素数：nyou
          nyou=nx*ny*nz
          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) yokoc
         write(7,500) tatec
         write(7,500) ziku
c
c 境界条件
         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 上下対称面（zx面）のy方向変位を拘束
        do k=1,nz
         write(7,1000) (nx+1)*(ny+1)*k+1
         write(7,1000) (nx+1)*(ny+1)*k+nx+1
         write(7,500) 1.d0
         write(7,500) 0.d0
         write(7,500) 1.d0
        end do
c 固定端の左右対称面(y軸)のx,z方向変位を固定
        do 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
        end do
c 左右対称面(yz面)のx方向変位を固定
       do  k=1,nz
        do 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
        end do
       end do
c 対称線（z軸）のxy方向変位を拘束
        do k=1,nz
         write(7,1000) (nx+1)*(ny+1)*k+1
         write(7,1000) (nx+1)*(ny+1)*k+1
         write(7,500) 0.d0
         write(7,500) 0.d0
         write(7,500) 1.d0
        end do
c 対称点（原点）を固定
         write(7,1000) 1
         write(7,1000) 1
         write(7,500) 0.d0
         write(7,500) 0.d0
         write(7,500) 0.d0
c
c
c *********** 載荷条件 ************
c
c 境界条件に合わせて温度荷重を与えるのは面倒なので
c 全節点に温度荷重を与えてしまって、mage2.f で
c 境界条件を入れる時に、拘束節点の荷重ベクトルの成分を
c 0にする。
c 全荷重ベクトルfそのものをtyokuon.fに渡すので、
c 載荷箇所とか載荷開始点終了点は要らない。
          do 20 i=0,nset-1
         write(7,500) f(i*3+1)
         write(7,500) f(i*3+2)
         write(7,500) f(i*3+3)
20          continue
c
        do i=1,3
        do j=1,3
         write(7,900) a(i,j)
        end do
        end do
         write(7,500) gxya
         write(7,500) gxza
         write(7,500) gyza
c
        do i=1,3
        do j=1,3
         write(7,900) b(i,j)
        end do
        end do
         write(7,500) gxyb
         write(7,500) gxzb
         write(7,500) gyzb
c
        do i=1,3
        do j=1,3
         write(7,900) c(i,j)
        end do
        end do
         write(7,500) gxyc
         write(7,500) gxzc
         write(7,500) gyzc
c
c
         write(7,500) df
         write(7,1000) nhikak
c 端部の図心変位のみを出す場合
         write(7,1000) ngamen
         write(7,1000) ngame1
c 端部の全節点の変位をファイル出力
         write(7,1000) nfairu
        do i=nfair1,nset
         write(7,1000) i
        end do
         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 yokoc*tateaの接着剤：2
c yokob*tateaの集成材：3
c yokoa*tatebの集成材：4
c yokoc*tatebの集成材：5
c yokob*tatebの集成材：6
c
         subroutine bunrui(nx,ny,nz,nxa,nya,nxc,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 左下の鋼板要素：横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(全部鋼材)だったら、ここで戻る
          if( (nxa.eq.nx).and.(nya.eq.ny) ) return
c
c 鋼板右隣の接着剤要素：横(nxc)*縦(nya)分割
           do i=nxa+1,nxa+nxc
              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
c 鋼板右隣の木材要素：横(nx-nxc-nxa)*縦(nya)分割
              iac=nxa+nxc
           do i=iac+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)=3
            end do
            end do
            end do
c もし鋼板が上から下までぶっとおっていたら、ここで戻る
         if(nya.eq.ny) return
c
c 鋼板上の左端の木材要素：横nx*縦(ny-nya)分割
           do i=1,nxa
           do j=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
c
c 接着剤要素の上の木材要素：横(nxc)*縦(ny-nya)分割
           do i=nxa+1,nxa+nxc
           do j=1,ny-nya
           do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=5
           end do
           end do
           end do
c
100        continue
c 残りの木材要素：横(nx-nxa-nxc)*縦(ny-nya)分割
              iac=nxa+nxc
           do i=iac+1,nx
           do j=1,ny-nya
           do k=1,nz
              ijk=nx*ny*(k-1)+nx*(j-1)+i
              ijkab(ijk)=6
           end do
           end do
           end do
          return
          end
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 更にbtをかけたもの(bt*e66*ep)を温度荷重にする
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,epx,epy,epz,fab,nab,bt)
         implicit real*8 (a-h, o-z)
         dimension fab(6,24),bt(24,6),e66(6,6),ep(6),ee(6)
          do i=1,24
             fab(nab,i)=0.d0
          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
c 温度ひずみ
          ep(1)=epx
          ep(2)=epy
          ep(3)=epz
          ep(4)=0.d0
          ep(5)=0.d0
          ep(6)=0.d0
c
         do i=1,6
            ee(i)=0.d0
            ee(i+3)=0.d0
         do j=1,6
            ee(i)=ee(i)+e66(i,j)*ep(j)
         end do
         end do
c
         do i=1,24
            fab(nab,i)=0.d0
         do j=1,6
            fab(nab,i)=fab(nab,i)+bt(i,j)*ee(j)
         end do
         end do
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
