! CalculiX用のinpファイルを作るためのプログラム（C3D4四面体要素用）
! まず、Salomeで適当に直方体要素のメッシュを切る
! それをUNV形式でsalome1.unvに保存。
! 拘束条件と載荷条件を適切に設定。
!ccxc3d4unv.f90をコンパイルして実行するとccxc3d4unv.inpが出来る。
!
implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000)

open(7,file='daiya2020.unv')
open(8,file='daiya2020.inp')
ell=1.0d-1
nsettensuu=34952  !節点数を入力
nedgeyouso=10600    !1D(edges)要素の数の入力
nfaceyouso=68968 !2D(faces)要素の数の入力
nvolumeyouso=106894 !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)
waru=1.d3
write(8,'(i6,a,1pe13.5,a,1pe13.5a,1pe13.5)') nset,',',x(i)/waru&
&,',',y(i)/waru,',',z(i)/waru
end do
write(8,'("** 要素番号とその要素を囲む節点番号")')
write(8,'("*ELEMENT, TYPE=C3D4, ELSET=Eall")')
read (7,'()') 
read (7,'()') 
read (7,'()') 
do i=1,nedgeyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
do i=1,nfaceyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
end do
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
end do
!read (7,'()') 
!read (7,'()') 
!read (7,'()') 
!write(8,'("**ここから節点グループをよんで飛ばす")')
!read (7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!!節点数が偶数の場合
!if (mod (nyou8,2) == 0)   then
!do i=1,nyou8/2+1
!read (7,'()') 
!end do
!!節点数が奇数の場合
!else
!do i=1,(nyou8+1)/2+1
!read (7,'()') 
!end do
!end if 
!read (7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!if (mod (nyou8,2) == 0)   then!節点数が偶数の場合
!do i=1,nyou8/2+1
!read (7,'()') 
!end do
!else!節点数が奇数の場合
!do i=1,(nyou8+1)/2+1
!read (7,'()') 
!end do
!end if 
!!write(8,'("**ここから1個目の要素グループをよんでかく")')
!write(8,'("** 要素グループを以下に挿入")')
!write(8,'("*ELSET,ELSET=E1")')
!read (7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!if (mod (nyou8,2) == 0)   then!要素数が偶数の場合
!read (7,'()') 
!do i=1,nyou8/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2,','
!write(8,*) n6,','
!end do
!else!要素数が奇数の場合
!read (7,'()') 
!do i=1,(nyou8-1)/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2,','
!write(8,*) n6,','
!end do
!read(7,*) n1,n2,n3,n4
!write(8,*) n2,','
!end if
!write(8,'("**ここから２個目の要素グループをよんでかく")')
!write(8,'("*ELSET,ELSET=E2")')
!read (7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!if (mod (nyou8,2) == 0)   then!要素数が偶数の場合
!read (7,'()') 
!do i=1,nyou8/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2,','
!write(8,*) n6,','
!end do
!else!要素数が奇数の場合
!read (7,'()') 
!do i=1,(nyou8-1)/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2,','
!write(8,*) n6,','
!end do
!read(7,*) n1,n2,n3,n4
!write(8,*) n2,','
!end if
close(7)
!write(8,'("** 要素グループここまで")')
!write(8,'("** 拘束節点を以下に挿入")')
write(8,'("*NSET,NSET=kousoku")')
do i=1,nsettensuu    !節点数
if ( -0.0000000000000001E+00<z(i).and.z(i)<0.0000000000000001E+00) then!拘束点いれる
write(8,*) i,','
end if
end do
write(8,'("** 拘束節点ここまで")')
write(8,'("** 載荷節点を以下に挿入")')
write(8,'("*NSET,NSET=saika")')
do i=1,nsettensuu    !節点数
if (102.0<z(i).and.z(i)<102.491997) then!載荷条件入れる102.491997
write(8,*) i,','
end if
end do
write(8,'("** 拘束節点ここまで")')
write(8,'("** 載荷条件を以下に挿入")')
write(8,'("** 載荷条件ここまで")')
write(8,'("** 材料の設定")')
write(8,'("*MATERIAL,NAME=EL")')
write(8,'("*ELASTIC")')
write(8,'("2.99d9,0.3")')!ヤング率
write(8,'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
write(8,'("** 拘束条件を以下に挿入")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku,1,3")')!拘束条件
write(8,'("** 拘束条件ここまで")')
write(8,'("** 載荷条件を以下に挿入")')
write(8,'("*CLOAD")')
write(8,'("saika,3,-0.005")')!saikaという条件でx=1,y=2,z=3
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
