<br />
<b>Deprecated</b>:  The each() function is deprecated. This message will be suppressed on further calls in <b>/home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php</b> on line <b>456</b><br />
! 2012/11/2版
! Salomeで生成したメッシュデータをGmshでinpに変換して
! CalculiX用のinpファイルを作るためのプログラム（C3D8直方体要素用）
! まず、Salomeで適当に直方体要素のメッシュを切る
! それをUNV形式でsalome.unvに保存。
! 次にGmshのFileのMergeでsalome.unvを読み込んで、
! Abaqus inp形式でsalome.inpに保存。
! salome.inpの冒頭にある節点番号と座標が並んでいる部分を
! setten.inpにコピーし、冒頭に節点数を記入。
! salome.inpの末尾にある要素番号と要素を囲む節点番号が並んでいる部分を
! youso.inpにコピーし、冒頭に要素数を記入。
! ccxs6を実行すると、ccxs6.inpが生成されているので、
! cgx -c ccxs6.inp で読み込んで拘束節点と載荷節点を設定し、それらを出力。
! 出力した拘束節点と載荷節点をccxs6.inpに読み込んで
! 拘束条件と載荷条件を適切に設定。
!
program ccxs6
implicit real*8(a-h,o-z)
dimension x(3000000),y(3000000),z(3000000),mesh(1000000,6),narabi(10000,10000)
!実行時にセグメンテーション違反が出る場合は、setten.inp内の節点数を確認し、
!x,y,zやnarabiの次元を増やすと解決する場合があるが、コンパイルが遅くなる。
!dimension x(3000000),y(3000000),z(3000000),mesh(1000000,6),narabi(16000,16000)
open(7,file='setten.inp')
read(7,*) nset
do i=1,nset
!read(7,*) nset,x,y,z
!write(8,*) nset,',',x,',',y,',',z
 read(7,*) kazu,x(i),y(i),z(i) !各節点の座標値の読み込み
end do
close(7)
!
open(7,file='youso.inp')
read(7,*) nyou
!
!
do i=1,nyou
!read(7,*) nyou,n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) i,',',n1,',',n2,',',n3,',',n4,',',n5,',',n6,',',n7,',',n8
read(7,*) kazu,mesh(i,1),mesh(i,2),mesh(i,3)
end do
close(7)
!
!
do i=1,nset
do j=1,nset
!隣り合う節点i,jの中間の節点番号の初期値をすべて-1にしておく
narabi(i,j)=-1
end do
end do
!
!
!中間点に新たに節点番号をふっていくのだが、
!頂点だけの節点数がnset個で、blenderの節点番号は0から始まるから、
!まだ使ってない一番若い節点番号はnset-1+1番
!2012/8/9メモ
! iset=nset-1 !こうしておけば、最初のiset+1=nsetとなる。
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要素目の頂点4に覚えさせる
  !頂点4の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
!
!
!!!!!!!!!!!!!!!!!!!!!!
print'("** blenderのobjファイルをccxのS6要素inp化")'
print'("**")'
print'("*NODE, NSET=Nall")'
do i=1,nset
   !print'(i6,a,f9.3,a,f9.3,a,f9.3)'&
   print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
   !print* &
   &,i,',',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),',',mesh(i,2),',',mesh(i,3),',',mesh(i,4)&
&   ,',',mesh(i,5),',',mesh(i,6)
end do
!拘束条件や材料条件はSTEP行よりも上に書かないとエラーが出る
!拘束条件を手動で書きこむ場合
!print'("*NSET,NSET=Nkousoku")'
!print'("** ここに拘束したい節点番号とコンマを1行ずつ並べる")'
!print'("** 例えば拘束節点が1と2だったら")'
!print'(i6,a)',1,','
!print'(i6,a)',2,','
!print'("** みたいに")'
print'("** 拘束条件をcgxで生成したkousoku.namから読み込む")'
write(8,'("** cgxで形状を見る場合は以下は**INCLUDEでコメント")')
write(8,'("**INCLUDE, INPUT=kousoku.nam")')
write(8,'("** ccxで計算するときは上の**INCLUDEの*を1個に")')
print'("*BOUNDARY")'
print'("** 以下のNkousoku,に続けて拘束変位1=x,2=y,3=zをコンマ区切りで")'
print'("Nkousoku,1,3")'
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC")'
e=210000. ! ヤング率(MPa)
poi=0.3   !  ポアソン比
print*, e,',',poi
t=10.d-3 !板の厚さ(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行ずつ並べる")'
!print'("** 例えば載荷節点が3と4だったら")'
!print'(i6,a)',3,','
!print'(i6,a)',4,','
!print'("** みたいに")'
print'("** 載荷条件をcgxで生成したsaika.namから読み込む")'
write(8,'("** cgxで形状を見る場合は以下は**INCLUDEでコメント")')
write(8,'("**INCLUDE, INPUT=saika.nam")')
write(8,'("** ccxで計算するときは上の**INCLUDEの*を1個に")')
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
p=1.d-0 !載荷荷重(MN)
nsaika=2 !載荷節点数
print'("** 以下のNsaikaに続けて載荷方向、載荷荷重を書く")'
print'("** 例えば載荷方向が3=zで、載荷荷重がp/nsaikaなら")'
print*,'Nsaika,3,',p/nsaika
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='ccxs6.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
!
!
close(8)
end