c 直方６面体立体要素の有限要素プログラム tyoku.f で
c z軸に横たわる片持ち梁の自由端に引張り（z方向）荷重を与える問題を
c 1/4 の角柱で解く（1/4解析）入力データを作るためのプログラム。
c このプログラムを実行させると、hipp.d というデータができあがるから
c tyoku<kata.d （hoge というファイルに出力するなら tyoku<kata.d>hoge）
c と実行させればよい。
c
c 02/8/1版
c プログラムを作った人：後藤文彦
c 改造などは自由ですが、もとのプログラムが、
c http://gthmhk.virtualave.net/programoj/tyoku/hipp.f
c であることをプログラム中に明記して戴けると嬉しいです。
c バグの報告は、
c http://gthmhk.virtualave.net/denbin.html
c にお願いします。
c
        program hipp
           implicit real*8 (a-h, o-z)
           dimension bc(9400*3), f(9400*3), a(3,3),
     & lkais(100),lryou(100)
c 梁の断面の横の長さ:yoko, 縦の長さ:tate
            yoko=36.d0
            tate=36.d0
c 梁の軸長：ziku
            ziku=36.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 節点数：nset
           nset=(nx+1)*(ny+1)*(nz+1)
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 載荷箇所の数：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 まず最初に、自由端の全ての節点に fz=2.d0 を与えてしまう。
c 1箇所目の載荷箇所における載荷開始節点番号:lkais(1)
c 自由端の節点番号の一番若い節点（左上のかど）
            lkais(1)=(nx+1)*(ny+1)*nz+1
c 1箇所目の載荷箇所における載荷終了節点番号:lryou(1)
            lryou(1)=nset
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=0.d0
             fz2=2.d0/fn
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=0.d0
              fz1=1.d0/fn
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=0.d0
             fz4=4.d0/fn
            endif
c
c ******* 材料定数 *********
c ヤング率：exx,eyy,ezz
            exx=16.d0/3.d0
            eyy=16.d0/3.d0
            ezz=16.d0/3.d0
c ポアソン比：poixy,poiyz,poizx
            poixy=1.d0/3.d0
            poixz=1.d0/3.d0
            poiyz=1.d0/3.d0
c せん断弾性係数：gxy,gyz,gzx
            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(2,1)
            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=1
c
c 変位を画面出力したい節点の数：ngamen
c 例えば、先端の節点すべてだったら：
            ngamen=(nx+1)*(ny+1)
c 1つ目の出力したい節点番号：ngame1
            ngame1=(nx+1)*(ny+1)*nz+1
c 変位をファイル出力したい節点の数：nfairu
            nfairu=1
c 1つ目の出力したい節点番号：nfair1
c 例えば、先端の中央の節点の節点番号だったら：
c            nfair1=(nx+1)*(ny+1)*nz
c     &            +(nx+1)*int(float(ny)/2.d0)
c     &            +int(float(nx)/2.d0)+1
c 先端の対称軸（４倍モデルの図心）の節点番号だったら、
            nfair1=(nx+1)*(ny+1)*nz
     &            +1
c
c
c tyoku.f 用の入力データを kata.d というファイルに出力する
c
         open(7,file='hipp.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 上下対称面（zx面）のy方向変位を拘束
        do 60 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
60      continue
c 固定端の左右対称面(y軸)のx,z方向変位を固定
        do 70 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
70      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 対称線（z軸）のxy方向変位を拘束
        do 80 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
80      continue
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
         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
         write(7,1000) ngamen
        do 30 i=ngame1, nset
         write(7,1000) i
30      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
