c mozi-koodo: EUC
c 直方６面体立体要素の有限要素プログラム tyoku.f で
c z軸に沿って横たわる片持ち梁の自由端に集中荷重（y方向）を与える問題を
c yz面を対称面とする半解析で解かせるための入力データを作るプログラム。
c このプログラムを実行させると、mage.d というデータができあがるから
c tyoku<mage.d （ hogeにファイル出力するなら tyoku<mage.d>hoge ）
c と実行させればよい。
c
c 02/8/1版
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c このプログラムは無保証です。                           c
c このプログラムは改造／再配布して構いません。           c
c 但し、私的利用に留まらない場合は、もとのプログラムが   c
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/      c
c から取ってきたものと分かるように明示して下さい。       c
c 改造／再配布後も上記の条件を保持して下さい。           c
c バグの報告は、                                         c
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg   c
c にお願いします。                                       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
        program mage
           implicit real*8 (a-h, o-z)
           dimension a(3,3),lkais(100),lryou(100)
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 梁の断面の横の長さ:yoko, 縦の長さ:tate
            yoko=1.d0
            tate=2.d0
c 梁の軸長：ziku
            ziku=5.d0
c 断面横方向の要素分割数：nx
c 断面縦方向の要素分割数：ny
c 軸方向の要素分割数：nz
          write(*,*) '要素分割数を入力して下さい'
          write(*,*) '（横,縦,軸長の順にコンマ区切りで）'
          write(*,*) '(縦+1)×(横+1)×(軸長+1)<9400'
          write(*,*) '例えば、6,10,100 や 6,12,99 や 20,20,20 など'
          read(*,*) nx,ny,nz
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/fn を与えてしまう。
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 次に、自由端の周を除くすべての節点に fy=4.d0/fn を与えなおす。
           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 ヤング率：exx,eyy,ezz
            exx=2.d9
            eyy=2.d9
            ezz=2.d9
c ポアソン比：poixy,poixz,poiyz
            poixy=1.d0/3.d0
            poixz=1.d0/3.d0
            poiyz=1.d0/3.d0
c せん断弾性係数：gxy,gxz,gyz
            gxy=ezz/2.d0/(1.d0+poixy)
            gxz=gxy
            gyz=gxy
c
c aに柔性（コンプライアス）行列を入れる
            a(1,1)=1.d0/exx
            a(1,2)=-poixy/exx
            a(1,3)=-poixz/exx
            a(2,1)=a(1,2)
            a(2,2)=1.d0/eyy
            a(2,3)=-poiyz/eyy
            a(3,1)=a(1,3)
            a(3,2)=a(2,3)
            a(3,3)=1.d0/ezz
c
c 荷重増分：df
            write(*,*) '荷重増分を入力して下さい'
            read(*,*) df
c 荷重増分ステップ数：nstep
            write(*,*) '荷重増分ステップ数を入力して下さい'
            read(*,*) nstep
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='mage.d')
         write(7,500) yoko
         write(7,500) tate
         write(7,500) ziku
         write(7,1000) nx
         write(7,1000) ny
         write(7,1000) nz
         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) gxy
         write(7,500) gxz
         write(7,500) gyz
         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
