!14/11/6 滝田さんのc3d4unv.f90を体積を求めるように修正
!14/5/12 Salome6.6用にSalome内でmm単位で書くと、inpファイルでm単位
!14/5/12 載荷と荷重を追加。
!c3d4unv.f90にSalomeの.unvの最初のいらない部分を飛ばすread文を追加。
!CalculiX用のinpファイルを作るためのプログラム（C3D4四面体要素用）
!まず、Salomeで適当に直方体要素のメッシュを切る
!それをunv保存。
!open(7,file='Mesh_1.unv')のファイル名を読み込むファイル名に
!open(8,file='Mesh_1.inp')を吐き出すファイル名(何でもいい)にする。
!meshをcomputeした時に出る画面の
!節点数とedges,faces,volumeを適切に入力。
!拘束条件と載荷条件を適切に設定。
!f90をコンパイルして実行するとinpファイルが出来る。
implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000) !節点座標
open(7,file='Mesh_1.unv')
open(8,file='mesh_1.inp')
p=1.d-3!荷重(MN)
ell=0.8d0 !inpファイル内での梁の長さ(載荷節点用)
nsettensuu=2349  !節点数を入力(Node)
nedgeyouso=422  !1D(edges)要素の数の入力
nfaceyouso=3264  !2D(faces)要素の数の入力
nvolumeyouso=8274 !3D(volume)要素の数の入力
write(8,'("** 節点番号と座標")')
write(8,'("*NODE, NSET=Nall")')
do i=1,19
read (7,'()') 
end do
!節点番号と節点座標を読み込む
do i=1,nsettensuu    !節点数
read (7,*) nset,n1,n2,n3
read(7,*) x(i),y(i),z(i)
write(8,'(i6,a,1pe13.5,a,1pe13.5a,1pe13.5)') nset,',',x(i)/1.0d3&
&,',',y(i)/1.0d3,',',z(i)/1.0d3
end do
!
write(8,'("** 要素番号とその要素を囲む節点番号")')
write(8,'("*ELEMENT, TYPE=C3D4, ELSET=Eall")')
read (7,'()') 
read (7,'()') 
read (7,'()') 
!まず、Edges情報を読む（使わない）
do i=1,nedgeyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
!次に Faces情報を読む（使わない）
do i=1,nfaceyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
end do
!次に Volumes情報を読む（これを使う)
taiseki=0.d0
do i=1,nvolumeyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) nset1,nset2,nset3,nset4
write(8,*) nyou-nfaceyouso-nedgeyouso,',',nset1,',',nset2,',',nset3,',',nset4
! 体積の計算
x1=x(nset1);y1=y(nset1);z1=z(nset1)
x2=x(nset2);y2=y(nset2);z2=z(nset2)
x3=x(nset3);y3=y(nset3);z3=z(nset3)
x4=x(nset4);y4=y(nset4);z4=z(nset4)
taiseki=taiseki+(&
& (x4-x1)*((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)) &
&+(y4-y1)*((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)) &
&+(z4-z1)*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)) &
&)/6
end do
print*,'体積は',taiseki
close(7)
!write(8,'("** 要素グループここまで")')
!write(8,'("** 拘束節点を以下に挿入")')
write(8,'("*NSET,NSET=kousoku1")')
!拘束節点の座標
yk1=0.
zk1=100.
yk2=0.
zk2=900.
do i=1,nsettensuu    !節点数
if(yk1-1.d-3<y(i).and.y(i)<yk1+1.d-3  &
&  .and. zk1-1.d-3<z(i).and.z(i)<zk1+1.d-3) then
write(8,*) i,','
end if
end do
write(8,'("*NSET,NSET=kousoku2")')
do i=1,nsettensuu    !節点数
if(yk2-1.d-3<y(i).and.y(i)<yk2+1.d-3  &
&  .and. zk2-1.d-3<z(i).and.z(i)<zk2+1.d-3) then
write(8,*) i,','
end if
end do
write(8,'("** 拘束節点ここまで")')
write(8,'("** 載荷節点を以下に挿入")')
write(8,'("*NSET,NSET=saika")')
nsaikasettensuu=0
!載荷節点の座標
ys=20.
zs=500.
do i=1,nsettensuu    !節点数
if(ys-1.d-3<y(i).and.y(i)<ys+1.d-3  &
&  .and. zs-1.d-3<z(i).and.z(i)<zs+1.d-3) then
write(8,*) i,','
nsaikasettensuu=nsaikasettensuu+1
end if
end do
write(8,'("** 拘束節点ここまで")')
write(8,'("** 載荷条件を以下に挿入")')
write(8,'("** 載荷条件ここまで")')
write(8,'("** 材料の設定")')
write(8,'("*MATERIAL,NAME=EL")')
write(8,'("*ELASTIC")')
write(8,'("206.0d3,0.3")')  !ヤング率(MPa), ポアソン比
write(8,'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
write(8,'("** 拘束条件を以下に挿入")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku1,1,3")')
write(8,'("kousoku2,1,2")')
write(8,'("** 拘束条件ここまで")')
write(8,'("** 載荷条件を以下に挿入")')
write(8,'("*CLOAD")')
! "saika"で定義した節点の2(y)方向にp/載荷節点数を載荷
write(8,'(A,a,i6,a,1pe13.5)')'saika',',',2,',',p/nsaikasettensuu
write(8,'("** 載荷条件ここまで")')
write(8,'("** 出力の設定")')
write(8,'("*NODE PRINT,NSET=saika")')
write(8,'("U")')
write(8,'("*EL PRINT,ELSET=Eall")')
write(8,'("S")')
write(8,'("*NODE FILE")')
write(8,'("U")')
write(8,'("*EL FILE,POSITION=AVERAGED AT NODES")')
write(8,'("S")')
write(8,'("*END STEP")')
close(8)
end
