!14/8/6 B32要素用 SalomeでLineでトラスを作り,
!全て選択し１つのWireにする.
!その後,作ったWireを,
!1Dでメッシュを切るが、number of segmentsで分割数指定して切る.
!*注意*number of segmentsは3以上の偶数にする##############
!拘束載荷と荷重は各自設定が必要.
!14/5/12 Salome6.6用にSalome内でmm単位で書くと、inpファイルでm単位
!c3d4unv.f90にSalomeの.unvの最初のいらない部分を飛ばすread文を追加。
!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),e66(6,6),b(3,3)
open(7,file='tora4.unv')
open(8,file='tora4.inp')
p=-1.0d-3!荷重(MN)
nsettensuu=80!節点数を入力(Node)
nedgeyouso=80!1D(edges)要素の数の入力
!#######木材の直行異方性の設定############
!木材ヤング率の入力
ezz=6.000d3
exx=ezz/25.0d0
eyy=ezz/25.0d0
!ポアソン比
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!#################################
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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=B32, ELSET=Eall")')
do i=1,3
read (7,'()') 
end do
!!要素数が偶数
!if (mod (nedgeyouso,2) == 0) then
do i=1,nedgeyouso/2
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) ntora1,ntora2
read(7,*) nyou2,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) ntora2,ntora3
write(8,*) i,',',ntora1,',',ntora2,',',ntora3
end do
!!奇数
!else
!do i=1,(nedgeyouso-1)/2
!read(7,*) nyou,n1,n2,n3,n4,n5
!read(7,*) n6,n7,n8
!read(7,*) ntora1,ntora2
!read(7,*) nyou2,n1,n2,n3,n4,n5
!read(7,*) n6,n7,n8
!read(7,*) ntora2,ntora3
!write(8,*) i,',',ntora1,',',ntora2,',',ntora3
!end do
!read(7,*) nyou2,n1,n2,n3,n4,n5
!read(7,*) n6,n7,n8
!read(7,*) ntora4,ntora5
!write(8,*) i+1,',',ntora3,',',ntora4,',',ntora5
!end if
close(7)
!write(8,'("** 要素グループここまで")')
!write(8,'("** 拘束節点を以下に挿入")')
!write(8,'("*NSET,NSET=kousoku")')
!do i=1,nsettensuu    !節点数
!if (-0.0001<z(i).and.z(i)<0.0001) then
!write(8,*) i,','
!end if
!end do
!write(8,'("** 拘束節点ここまで")')
!write(8,'("** 載荷節点を以下に挿入")')
!write(8,'("*NSET,NSET=saika")')
!nsaikasettensuu=0
!do i=1,nsettensuu    !節点数
!if (ell-0.00002<z(i)/1.0d3.and.z(i)/1.0d3<ell+0.00002) then!片持ち梁の端っこのz座標を条件にする
!write(8,*) i,','
!nsaikasettensuu=nsaikasettensuu+1
!end if
!end do
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,'("*MATERIAL,NAME=EL")')
write(8,'("*ELASTIC,TYPE=ORTHOROPIC")')
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
    & gxz,gxy,poiyx,poizx,poizy)
write(8,'(f12.4,a,f12.4,a,f12.4,a,&
    &f12.4,a,f12.4,a,f12.4,a,f12.4,&
    & a ,f12.4,a)')&
    &e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
    &,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
write(8,'(f12.4,a,a)')e66(4,4),',','293'
!write(8,'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")')
write(8,'("*BEAM SECTION,ELSET=Eall,MATERIAL=EL ,SECTION=CIRC")')
write(8,'("0.1,0.1")')
write(8,'("0.,1.,0.")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
write(8,'("** 拘束条件を以下に挿入")')
write(8,'("*BOUNDARY")')
write(8,'("2,1,3")')
write(8,'("** 拘束条件ここまで")')
write(8,'("** 載荷条件を以下に挿入")')
write(8,'("*CLOAD")')
write(8,'(A,a,i6,a,1pe13.5)')'1',',',2,',',p
write(8,'("** 載荷条件ここまで")')
write(8,'("** 出力の設定")')
write(8,'("*NODE PRINT,NSET=Nall")')
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
!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
!  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz
!
!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
    &+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
    &-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!
!
return
end

