!mokuc3d10.f90をオンサイト木橋の座屈用にいらない部分を消してすっきりさせた。

! 拘束とか載荷の節点群は、座標指定ではなく、Salomeのmeshのグループで指定する。
! なので、Salomeで、Mesh_1を右クリック-->Create Groups from Geometry-->Nodes
! から、geometryで予めcreate groupしたグループを選択して作っておく。
! 材料別の要素群も同様に、meshのグループを作る必要があるが、
! 要素は並び順情報が要るので、Create Groups from Geometry-->Elements で選択する。

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

implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000)
dimension e66(6,6),b(3,3),mesh(30000000,10),ngroup(30000000)
dimension t(3,3),alp(1),bt(1),gmm(1)

!---各種入力データ---
!#--入出力ファイル名
open(7,file='unv/tamesi.unv')
open(8,file='tamesi/tamesi_yuki.inp')
!#--荷重
p=1.d-0/4!集中荷重/4(MN)
q=3.5d-3!積雪荷重(MN/m^3)
!#--寸法
ell=8395.d-3!橋長/2(m)
ell1=150.d-3!端から支点までの距離(m)
h=800.d-3!高さ(m)
b1=480.d-3!木材の幅員/2(m)
t1=9.d-3+100.d-3!鋼板の板厚+補剛材のフランジ幅(m)
!#--Mesh_1 の Mesh Information を見て入力
nnode=489610   !2次要素の節点数(Nodes)
nedge=9516     !1D(edges)の数
nface=118796    !2D(faces)の数
nvolume=280843 !3D(volumes)の数
!---sub mesh の Mesh Information を見て入力
!#--各部材の要素数(volumes)と固定軸回りのオイラー角？
nbuzai=1
pi=asin(1.d0)*2.d0
nyouso1=109821     !木材
alp(1)= 0.d0         ; bt(1)= 0.d0       ; gmm(1)=0.d0
nyouso0=171022    !鋼材
!#--材料定数
!---鋼材か樹脂のヤング率
young=206.d3
!young=2.84d3
!---鋼材か樹脂のポアソン比
poisson=0.3d0
!---木材のヤング率
ezz=7.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
!---

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

!---inpファイルの要素情報を読み取る---
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情報の読み飛ばし
do i=1,nface
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8,n9,n10,n11
end do
!---Volumes情報の読み込みと書き出し、ついでに体積計算
taiseki=0.d0
do i=1,nvolume
read(7,*) ngroup(i),n1,n2,n3,n4,n5
read(7,*) mesh(ngroup(i),1),mesh(ngroup(i),5),mesh(ngroup(i),2),mesh(ngroup(i),6)&
    &,mesh(ngroup(i),3),mesh(ngroup(i),7),mesh(ngroup(i),8),mesh(ngroup(i),9)
read(7,*) mesh(ngroup(i),10),mesh(ngroup(i),4)
!---体積の計算
x1=x(mesh(ngroup(i),1));y1=y(mesh(ngroup(i),1));z1=z(mesh(ngroup(i),1))
x2=x(mesh(ngroup(i),2));y2=y(mesh(ngroup(i),2));z2=z(mesh(ngroup(i),2))
x3=x(mesh(ngroup(i),3));y3=y(mesh(ngroup(i),3));z3=z(mesh(ngroup(i),3))
x4=x(mesh(ngroup(i),4));y4=y(mesh(ngroup(i),4));z4=z(mesh(ngroup(i),4))
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
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---

!---inpファイルの節点群(*NSET)を作る---
!---幅方向(x)対称面
write(8,'("*NSET,NSET=koteix")')
do i=1,nnode
if (x(i)/1.d3<b1+1.d-5.and.x(i)/1.d3>b1-1.d-5) then
write(8,*) i,','
endif
enddo
!---支点(y方向固定)
write(8,'("*NSET,NSET=koteiy")')
do i=1,nnode
if (y(i)/1.d3<h+1.d-5.and.y(i)/1.d3>h-1.d-5.and.&
   &z(i)/1.d3<ell1+1.d-5.and.z(i)/1.d3>ell1-1.d-5) then
write(8,*) i,','
endif
enddo
!---幅方向(z)対称面
write(8,'("*NSET,NSET=koteiz")')
do i=1,nnode
if (z(i)/1.d3<ell+1.d-5.and.z(i)/1.d3>ell-1.d-5) then
write(8,*) i,','
endif
enddo
!---梁中央の上縁
nsensaika=0
write(8,'("*NSET,NSET=sensaika")')
do i=1,nnode
if (z(i)/1.d3<ell+1.d-5.and.z(i)/1.d3>ell-1.d-5.and.&
   &y(i)/1.d3<1.d-5.and.y(i)/1.d3>-1.d-5) then
nsensaika=nsensaika+1
write(8,*) i,','
endif
enddo
!---梁の上面
nmensaika=0
write(8,'("*NSET,NSET=mensaika")')
do i=1,nnode
if (y(i)/1.d3<1.d-5.and.y(i)/1.d3>-1.d-5) then !上面に面載荷
nmensaika=nmensaika+1
write(8,*) i,','
endif
enddo

!---inpファイルの要素群(*ELSET)を作る---
nyou=0
!   以下のグループの並び順は、Salomeで sub mesh を作った順に並び替える
!#--木材
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E1")')
do i=1,int(nyouso1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso1,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!#--鋼材
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E0")')
do i=1,int(nyouso0/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso0,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!---
close(7)

!#--境界条件(拘束方向は、x:1,y:2,z:3)
write(8,'("*BOUNDARY")')
write(8,'("koteix,1")')
write(8,'("koteiy,2")')
write(8,'("koteiz,3")')

!#--材料の設定0
write(8,'("*MATERIAL,NAME=Ekou")')
write(8,'("*ELASTIC")')
write(8,'(1pe17.9,a,1pe17.9)')young,',',poisson
write(8,'("*SOLID SECTION,ELSET=E0,MATERIAL=Ekou")')
!   ここから、異方性の設定
!#--材料の設定:木材1
nb=1!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*MATERIAL,NAME=Emoku")')
write(8,'("*ELASTIC,TYPE=ORTHO")')
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)'),e66(4,4)
write(8,'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E1,MATERIAL=Emoku,ORIENTATION=OR1")')

!---計算の設定(*STEP~*END STEP)
write(8,'("*STEP")')
write(8,'("*STATIC")')
!#--荷重条件(載荷方向は、x:1,y:2,z:3)
write(8,'("*BUCKLE")')
write(8,'("10,0.01")')!出力モード数、精度
write(8,'("*CLOAD")')
!write(8,'(A,a,i6,a,1pe25.13)')'sensaika',',',2,',',p/nsensaika!3点曲げの場合
write(8,'(A,a,i6,a,1pe25.13)')'mensaika',',',2,',',q*(b1+t1)*ell/nmensaika!積雪荷重の場合
!#--出力の設定
write(8,'("*NODE PRINT,NSET=sensaika")')!梁中央の上縁の変位をたわみの確認用に出力する
write(8,'("U")')
!write(8,'("*NODE PRINT,NSET=Nall")')
!write(8,'("U")')
!write(8,'("*EL PRINT,GLOBAL=YES,ELSET=Eall")')
!write(8,'("S")')
write(8,'("*NODE FILE")')
write(8,'("U")')
write(8,'("*EL FILE,GLOBAL=YES")')
write(8,'("S,E")')
write(8,'("*END STEP")')
!---
close(8)
end

!-------- 直交異方性の定数を出すプログラム --------
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3),t(3,3),alp(1),bt(1),gmm(1)
!応力－ひずみ行列：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
!---固定軸回りのオイラー角？の回転行列
t=0.d0
t(1,1)=              cos(bt(nb))*cos(gmm(nb))
t(2,1)= sin(alp(nb))*sin(bt(nb))*cos(gmm(nb))+cos(alp(nb))*sin(gmm(nb))
t(3,1)=-cos(alp(nb))*sin(bt(nb))*cos(gmm(nb))+sin(alp(nb))*sin(gmm(nb))
t(1,2)=-cos(bt(nb))*sin(gmm(nb))                               
t(2,2)=-sin(alp(nb))*sin(bt(nb))*sin(gmm(nb))+cos(alp(nb))*cos(gmm(nb))
t(3,2)= cos(alp(nb))*sin(bt(nb))*sin(gmm(nb))+sin(alp(nb))*cos(gmm(nb))
t(1,3)= sin(bt(nb))                                                
t(2,3)=-sin(alp(nb))*cos(bt(nb))                                    
t(3,3)= cos(alp(nb))*cos(bt(nb))

return
end
