!Salomeで作った三角形板要素(2次要素)のunvファイルをCalculiX用inpファイルに変換

! 拘束とか載荷の節点群は、座標指定が楽かと思い、そうする。
! つまり、等価節点力にはなっていないので注意。

! ちなみに、  /#  で 修正必要箇所 に跳べる。

implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000) !節点座標

!---各種入力データ---
!#--ファイル名
open(7,file='s6.unv')
open(8,file='s6.inp')
!#--荷重
p=1.d-4!荷重(MN)
t=1.d-3!板厚(mm)
!#--Mesh_1 の Mesh Information を見て入力
nsetten=497   !2次要素の節点数(Nodes)
nedge=40    !1D(edges)の数
nface=228   !2D(faces)の数
nvolume=412 !3D(volume)の数
!---

!---inpファイルの節点情報(*NODE)を作る---
write(8,'("*NODE, NSET=Nall")')
do i=1,19
read (7,'()')!ヘッダーの読み飛ばし
end do
!---節点番号と節点座標の読み込みと書き出し
do i=1,nsetten
read (7,*) nset,n1,n2,n3
read(7,*) x(i),y(i),z(i)
write(8,'(i6,a,1pe17.9,a,1pe17.9a,1pe17.9)') nset,',',x(i)/1.0d3&
&,',',y(i)/1.0d3,',',z(i)/1.0d3
end do
!---

!---inpファイルの要素情報(*ELEMENT)を作る---
write(8,'("*ELEMENT, TYPE=S6, ELSET=Eall")')!今のところ1種材料
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---Edges情報の読み飛ばし
do i=1,nedge
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
!---Faces情報の読みこみと書き出し，ついでに面積計算
sall=0.d0
do i=1,nface
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) nset1,nset4,nset2,nset5,nset3,nset6
write(8,*) i,',',nset1,',',nset2,',',nset3,',',nset4,',',nset5,',',nset6
!---面積の計算
vx1=x(nset1)-x(nset3)
vx2=x(nset2)-x(nset3)
vy1=y(nset1)-y(nset3)
vy2=y(nset2)-y(nset3)
vz1=z(nset1)-z(nset3)
vz2=z(nset2)-z(nset3)
v1=vy1*vz2-vz1*vy2
v2=vz1*vx2-vx1*vz2
v3=vx1*vy2-vy1*vx2
s=sqrt(v1**2.d0+v2**2.d0+v3**2.d0)/2.d0
sall=sall+s
end do
print*,'面積は',sall
!---

!---inpファイルの節点群(*NSET)を作る---
!#--固定するグループ1の節点番号の読み込みと書き出し
write(8,'("*NSET,NSET=kotei1")')
do i=1,nsetten
if (x(i)/1.d3>-50.d-3-1.d-5.and.x(i)/1.d3<-50.d-3+1.d-5) then!x座標の範囲
write(8,*) i,','
endif
enddo
!#--固定するグループ2の節点番号の読み込みと書き出し
!write(8,'("*NSET,NSET=kotei2")')
!do i=1,nsetten
!if (x(i)/1.d3>-50.d-3-1.d-5.and.x(i)/1.d3<-50.d-3+1.d-5.and.&!x座標の範囲
!&y(i)/1.d3>-50.d-3-1.d-5.and.y(i)/1.d3<-50.d-3+1.d-5.and.&!y座標の範囲
!&z(i)/1.d3>-50.d-3-1.d-5.and.z(i)/1.d3<-50.d-3+1.d-5) then!z座標の範囲
!write(8,*) i,','
!endif
!enddo
!#--載荷するグループの節点番号の読み込みと書き出し
write(8,'("*NSET,NSET=saika")')
nsaika=0
do i=1,nsetten
if (x(i)/1.d3>50.d-3-1.d-5.and.x(i)/1.d3<50.d-3+1.d-5) then!x座標の範囲
write(8,*) i,','
nsaika=nsaika+1
endif
enddo
!---
close(7)

!---inpファイルに色々な情報を書く---
!#--材料の設定
write(8,'("*MATERIAL,NAME=EL")')
write(8,'("*ELASTIC")')
write(8,'("6.0d3,0.3")')!ヤング率(MPa), ポアソン比
write(8,'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")')
write(8,'(1pe13.5,A,A)')t,',,','EL'
!---計算の設定(*STEP~*END STEP)
write(8,'("*STEP")')
write(8,'("*STATIC")')
!#--拘束条件(拘束方向は、x:1,y:2,z:3,th_x:4,th_y:5,th_z:6)
write(8,'("*BOUNDARY")')
write(8,'("kotei1,1,6")')
!write(8,'("kotei2,1,3")')
!#--荷重条件(載荷方向は、x:1,y:2,z:3)
write(8,'("*CLOAD")')
write(8,'(A,a,i6,a,1pe17.9)')'saika',',',3,',',-p/nsaika
!---出力の設定
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")')
write(8,'("S")')
write(8,'("*END STEP")')
!---
close(8)
end
