! EUC-jp
! 文字コードはEUC-jp
! 自由形式のfortranで書いています。
! g77でコンパイルする場合は、
!g77 -ffree-form -o smv smv.f
! のようにコンパイルしてください。
! gfortran でコンパイルする場合は文字コードをutf-8にしてください。
!ble6ccx.fのプログラムを折り畳みの拘束、節点番号を出力させる
!ように書き換えたプログラム
!折り畳み円筒の座屈荷重を求めるのに必要な
!inpファイルを出力させるプログラム
! 09/3/22版
! プログラムを主に作った人：後藤文彦
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! このプログラムは無保証です。                           
! このプログラムは改造／再配布して構いません。           
! 但し、私的利用に留まらない場合は、もとのプログラムが   
! http://www.str.ce.akita-u.ac.jp/~kudo/prog/
! から取ってきたものと分かるように明示して下さい。       
! 改造／再配布後も上記の条件を保持して下さい。           
! バグの報告は、                                         
! http://www.str.ce.akita-u.ac.jp/~gotou/png/meeru.png
! にお願いします。                                       
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
program ble6ccx
implicit real*8 (a-h, o-z)
dimension x(3000000),y(3000000),z(3000000),&
&         mesh(1000000,6),narabi(10000,10000)
character mozi*20,mozi0*20
read*, mozi0 !1行目の文字を読み込む
read*, nset !節点数の読み込み
do i=0,nset-1
 read*, x(i),y(i),z(i) !各節点の座標値の読み込み
end do
!
open(7,file='yousosuu.txt') !yousosuu.fで数えた要素数ファイル
read(7,*) nyou !要素数の読み込み
read(7,*) nsyuubun,ntakabun !周方向、高さ方向の分割数の読み込み
do i=1,nyou
!各要素の頂点数?、各頂点の節点番号、カラーコード?の読み込み
read*,kazu,mesh(i,1),mesh(i,2),mesh(i,3),mozi
end do
!
do i=0,nset
do j=0,nset
!隣り合う節点i,jの中間の節点番号の初期値をすべて-1にしておく
narabi(i,j)=-1
end do
end do
!
!
!中間点に新たに節点番号をふっていくのだが、
!頂点だけの節点数がnset個で、blenderの節点番号は0から始まるから、
!まだ使ってない一番若い節点番号はnset-1番
iset=nset-1
do i=1,nyou !1要素目からnyou要素まで順番に中間節点をふる
 !頂点1-2間に中間節点番号がまだふられていなければ、
 if(narabi(mesh(i,1),mesh(i,2))==-1) then
  iset=iset+1 !節点番号を1増やして
  narabi(mesh(i,1),mesh(i,2))=iset !頂点1-2間と
  narabi(mesh(i,2),mesh(i,1))=iset !頂点2-1間とに節点番号isetをふり
  mesh(i,4)=iset !その節点番号をi要素目の頂点5に覚えさせる
  !頂点5のxyz座標を頂点1-2間の中点として与える
  x(iset)=(x(mesh(i,1))+x(mesh(i,2)))/2.d0
  y(iset)=(y(mesh(i,1))+y(mesh(i,2)))/2.d0
  z(iset)=(z(mesh(i,1))+z(mesh(i,2)))/2.d0
 else
  !頂点1-2間にすでに節点番号がふられていれば、
  mesh(i,4)=narabi(mesh(i,1),mesh(i,2)) !それを頂点5の節点番号にする
 end if
!       頂点2-3,3-1間も以下同様に
 if(narabi(mesh(i,2),mesh(i,3))==-1) then
  iset=iset+1
  narabi(mesh(i,2),mesh(i,3))=iset
  narabi(mesh(i,3),mesh(i,2))=iset
  mesh(i,5)=iset
  x(iset)=(x(mesh(i,2))+x(mesh(i,3)))/2.d0
  y(iset)=(y(mesh(i,2))+y(mesh(i,3)))/2.d0
  z(iset)=(z(mesh(i,2))+z(mesh(i,3)))/2.d0
 else
  mesh(i,5)=narabi(mesh(i,2),mesh(i,3))
 end if
!
 if(narabi(mesh(i,3),mesh(i,1))==-1) then
  iset=iset+1
  narabi(mesh(i,3),mesh(i,1))=iset
  narabi(mesh(i,1),mesh(i,3))=iset
  mesh(i,6)=iset
  x(iset)=(x(mesh(i,3))+x(mesh(i,1)))/2.d0
  y(iset)=(y(mesh(i,3))+y(mesh(i,1)))/2.d0
  z(iset)=(z(mesh(i,3))+z(mesh(i,1)))/2.d0
 else
  mesh(i,6)=narabi(mesh(i,3),mesh(i,1))
 end if
 nset=iset
end do
!
!
!!!!!!!!!!!!!!!!!!!!!!
!open(9,file='ble6ccx.inp')
print'("** blenderのobjファイルをccxのS6要素inp化")'
print'("**")'
print'("*NODE, NSET=Nall")'
do i=0,nset
   print'(i6,a,f9.6,a,f9.6,a,f9.6)'&
   &,i+1,',',x(i),',',y(i),',',z(i)
end do
print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
do i=1,nyou
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
   i,',',mesh(i,1)+1,',',mesh(i,2)+1,',',mesh(i,3)+1,',',mesh(i,4)+1&
&   ,',',mesh(i,5)+1,',',mesh(i,6)+1
end do
print'("*NSET,NSET=FIX")'
print'("** ここに拘束したい節点番号とコンマを1行ずつ並べる")'
!
!ここで下端開口部の拘束する節点番号を入力
do i=0,nsyuubun-2,2 !周方向に
print*, i*3+1,','
print*, i*3+2,','
print*, i*3+6,','
end do
print*,nsyuubun*3/2*(ntakabun+1)+1,','
do i=1,nsyuubun*3/2-1
print*,nsyuubun*3/2*(ntakabun+1)+2+(4*i),','
end do

!
!ここで上端開口部の拘束する節点番号を入力
print'("*NSET,NSET=FIX2")'
print'("** ここに拘束したい節点番号とコンマを1行ずつ並べる")'
do i=1,nsyuubun*3/2
print*, nsyuubun*3*ntakabun/2+i,','
end do
do i=0,nsyuubun*3/2-2
print*, (nsyuubun*ntakabun*6)-nsyuubun/2-(nsyuubun-4)+(3*i),','
end do
print*, nsyuubun*3*(2*ntakabun+1),','
!
!ここで境界条件を入力
print'("*BOUNDARY")'
print'("** 以下のFIX,に続けて拘束変位1=x,2=y,3=zをコンマ区切りで")'
print'("FIX,1,3")'
print'("FIX2,1,2")'
!if(nsyuubun % 4==0) then
!print*, nsyuubun*2-(nsyuubun/2-1),',2,3'!yz拘束
!print'("63,2,3")'!yz拘束
!else if(nsyuubun==10) then
!print'("105,2,3")'!yz拘束
!else if(nsyuubun==14) then
!print'("147,2,3")'!yz拘束
!else
!print*, nsyuubun*10+nsyuubun/2,',2,3'!yz拘束
!end if
!
!ここで材料定数を定義
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC")'
e=1300. ! ヤング率(MPa)
poi=0.3   !  ポアソン比
print*, e,',',poi
t=4.0d-4 !板の厚さ(m)
print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL,OFFSET=0.")'
print*,t
!温度設定
!print'("*INITIAL CONDITIONS,TYPE=TEMPERATURE")'
!print'("Nall,273.")'
!節点ごとに厚さを変えたい時
!print'("*NSET,NSET=ATUSA")'
!do i=1,nset
!print'(i6,a)',i,','
!end do
!print'("*NODAL THICKNESS")'
!print*,'ATUSA,',t
!
!ここで載荷する節点番号を入力
print'("*NSET,NSET=SAIKA")'
print'("** ここに載荷したい節点番号とコンマを1行ずつ並べる")'
do i=1,nsyuubun*3/2
print*, nsyuubun*3*ntakabun/2+i,','
end do
do i=0,nsyuubun*3/2-2
print*, (nsyuubun*ntakabun*6)-nsyuubun/2-(nsyuubun-4)+(3*i),','
end do
print*, nsyuubun*3*(2*ntakabun+1),','
!
print'("*STEP")'
!print'("*STEP,NLGEON")'
print'("*STATIC,SOLVER=SPOOLES")'
!座屈荷重を求める場合は、以下のBUCKLEを入力
print'("*BUCKLE")'
print'("3,0.01")'
print'("*CLOAD")'
p=4.9E-6 !載荷荷重(MN)
print'("** 以下のSAIKAに続けて載荷方向、載荷荷重を書く")'
print'("** 例えば載荷方向が3=zで、載荷荷重がp/nsaikaなら")'
!荷重は載荷節点数で割る
print*,'SAIKA,3,',-(p/2/(nsyuubun*3/2))
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
print'("*EL PRINT,ELSET=Eall")'
print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'
!
!
!!!!!!!!!!!!!!!!!!!!!!
!
open(8,file='ble6ccx.obj')
write(8,'("3DG1")')
write(8,*) nset+1
do i=0,nset
write(8,*) x(i),y(i),z(i)
end do
!
do i=1,nyou
   write(8,*) 6,mesh(i,1),mesh(i,2),mesh(i,3),mesh(i,4)&
&              ,mesh(i,5),mesh(i,6),' ','0xcccccc'
end do
!
end

