c 直方６面体立体要素の有限要素プログラム tyoku.f で
c z軸に横たわる片持ち梁の自由端に引張り（z方向）荷重を与える問題を
c 1/4 の角柱で解く（1/4解析）入力データを作るためのプログラム。
c このプログラムを実行させると、hipp.d というデータができあがるから
c tyokuhoge）
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