! 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