!gom.f90を改良
!節点座標の出力桁を多くして、
!ゴムを小さくした時に節点座標が出力桁以下になり、
!四捨五入されてしまいメッシュの大きさが
!均一にならないことがある不具合を修正。
implicit real*8(a-h,o-z)
dimension x(315500),y(315500),z(315500),e66(6,6),b(3,3)&
&,meshmoku(300000,8),meshkou(300000,8),meshgom(300000,8)
!材料の寸法の設定
b1=1.2d-1!角材の幅
b2=9.0d-3!鋼板の幅
h1=1.2d-1!角材の高さ
h2=5.0d-1!鋼板の高さ
h3=h2-h1*2!間の鋼板の高さ
ell=7.0d0/2.0!スパンの半分の長さ
p=5.000d-2
! x,y,z方向の要素分割数
nx1=6!木材の分割数
nx2=6!鋼板の分割数
ny1=6!上と下部分の高さ方向の分割数
ny2=4!間の鋼板の高さ方向の分割数
nz=20!軸方向の分割数
!材料の数
nk=7!角材の数
nmk=2!鋼板の数
! 1要素のx,y,z方向の長さ：xn,yn,zn
xn1=b1/nx1!木材の1要素あたりの幅方向長さ
xn2=b2/nx2!鋼板の1要素あたりの幅方向長さ
yn1=h1/ny1!上と下部分の1要素あたりの高さ方向の長さ 
yn2=h3/ny2!間の鋼板の1要素あたりの高さ方向の長さ
zn=ell/nz!軸方向の1要素あたりの長さ
!ゴム系
b3=1.0d-5!ゴムの幅
nx3=4!ゴムの分割数
xn3=b3/nx3
!ヤング率の入力
ezz=6.921d3!z軸方向のヤング率
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!分割数
!一面とかよく使う記述を定義する
nx=(2*nx2+nx3*(nk+1)+(nx1)*nk+1)
nitimen=nx*2*(ny1+1)+(nx2+1)*(ny2-1)*2
! 下と上部分の節点　鋼板
do k=0,nz
do m=0,1
do l=0,1
do j=0,ny1
do i=1,nx2+1
nset=i+j*nx&
    &+l*(nx2+nx3*(nk+1)+nx1*nk)&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
x(nset)=xn2*(i-1)+l*(b2+b1*nk+b3*(nk+1))
y(nset)=yn1*j+m*(h1+h3)
z(nset)=zn*k
!print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
end do
! 鋼板 中部分の節点
do k=0,nz
do j=0,ny2-2
do l=0,1
do i=1,nx2+1
nset=nx*(ny1+1)+i+l*(nx2+1)+j*(2*(nx2+1))&
    &+k*nitimen
x(nset)=xn2*(i-1)+l*(b2+b1*nk+b3*(nk+1))
y(nset)=yn2*j+h1+yn2
z(nset)=zn*k
!print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
! ゴム
do k=0,nz
do m=0,1
do j=0,ny1
do l=0,nk
do i=1,nx3-1
nset=nx2+1+i+l*(nx3+nx1)+j*nx&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
x(nset)=xn3*(i-1)+b2+xn3+l*(b1+b3)
y(nset)=yn1*j+m*(h1+h3)
z(nset)=zn*k
!print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
end do
! 角材
do k=0,nz
do m=0,1
do j=0,ny1
do l=0,nk-1
do i=1,nx1+1
nset=nx2+nx3+i+l*(nx3+nx1)+j*nx&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
x(nset)=xn1*(i-1)+b2+b3+l*(b1+b3)
y(nset)=yn1*j+m*(h1+h3)
z(nset)=zn*k
!print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
end do
!i5~6が整数,aが文字,1pe17.9が小数.1peとかはEの＊乗
!接点番号と座標の出力
print'("*NODE, NSET=Nall")'
do i=1,(nz+1)*nitimen
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
&i,',',x(i),',',y(i),',',z(i)
nsettensuu=i
end do
!print*,nx
!print*,nitimen
!youso
!要素番号とその要素の8接点の接点番号
!上部分と下部分の鋼板
nyou=0
do k=0,nz-1 
do m=0,1
do l=0,1
do j=0,ny1-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=i+nx*j&
    &+l*(nx2+nx3*(nk+1)+nx1*nk)&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
meshkou(nyou,2)=meshkou(nyou,1)+nx
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
end do
end do
!真ん中の鋼板(最初の一段目)
do k=0,nz-1
do l=0,1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=nx*ny1+i+l*(nx-(nx2+1))&
    &+k*nitimen
meshkou(nyou,2)=nx*(ny1+1)+i+l*(nx2+1)&
    &+k*nitimen
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
!真ん中の鋼板(2段目から最後の前まで)
do k=0,nz-1
do j=0,ny2-3!ny2>3
do l=0,1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=nx*(ny1+1)+i&
    &+l*(nx2+1)+j*2*(nx2+1)+k*nitimen
meshkou(nyou,2)=meshkou(nyou,1)+(nx2+1)*2
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
end do
!真ん中の鋼板(最後の段)
do k=0,nz-1
do l=0,1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=nx*(ny1+1)+(ny2-2)*2*(nx2+1)&
    &+i+l*(nx2+1)&
    &+k*nitimen
meshkou(nyou,2)=nx*(ny1+1)+(ny2-1)*2*(nx2+1)&
    &+i&
    &+l*(nx-(nx2+1))&
    &+k*nitimen
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
!下と上部分 木材
nyou=0
do k=0,nz-1
do m=0,1
do j=0,ny1-1
do l=0,nk-1
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=nx2+nx3+i+l*(nx1+nx3)&
    &+j*nx&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
meshmoku(nyou,2)=meshmoku(nyou,1)+nx
meshmoku(nyou,3)=meshmoku(nyou,2)+nitimen
meshmoku(nyou,4)=meshmoku(nyou,1)+nitimen
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,2)+1
meshmoku(nyou,7)=meshmoku(nyou,3)+1
meshmoku(nyou,8)=meshmoku(nyou,4)+1
end do
end do
end do
end do
end do
!下と上部分 gom
nyou=0
do k=0,nz-1
do m=0,1
do j=0,ny1-1
do l=0,nk
do i=1,nx3
nyou=nyou+1
meshgom(nyou,1)=nx2+i+l*(nx3+nx1)&
    &+j*nx&
    &+m*((ny1+1)*nx+(nx2+1)*(ny2-1)*2)&
    &+k*nitimen
meshgom(nyou,2)=meshgom(nyou,1)+nx
meshgom(nyou,3)=meshgom(nyou,2)+nitimen
meshgom(nyou,4)=meshgom(nyou,1)+nitimen
meshgom(nyou,5)=meshgom(nyou,1)+1
meshgom(nyou,6)=meshgom(nyou,2)+1
meshgom(nyou,7)=meshgom(nyou,3)+1
meshgom(nyou,8)=meshgom(nyou,4)+1
end do
end do
end do
end do
end do
print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")'
do i=1,2*nx2*(ny1*2+ny2)*nz
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i,',',&
&meshkou(i,1),',',&
&meshkou(i,2),',',&
&meshkou(i,3),',',&
&meshkou(i,4),',',&
&meshkou(i,5),',',&
&meshkou(i,6),',',&
&meshkou(i,7),',',&
&meshkou(i,8)
end do
nyousosuukou=2*nx2*(ny1*2+ny2)*nz
print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")'
do i=1,2*nx1*ny1*nk*nz
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyousosuukou,',',&
&meshmoku(i,1),',',&
&meshmoku(i,2),',',&
&meshmoku(i,3),',',&
&meshmoku(i,4),',',&
&meshmoku(i,5),',',&
&meshmoku(i,6),',',&
&meshmoku(i,7),',',&
&meshmoku(i,8)
end do
nyousosuumoku=2*nx1*ny1*nk*nz
print'("*ELEMENT, TYPE=C3D8, ELSET=Egom")'
do i=1,2*nx3*ny1*(nk+1)*nz
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyousosuukou+nyousosuumoku,',',&
&meshgom(i,1),',',&
&meshgom(i,2),',',&
&meshgom(i,3),',',&
&meshgom(i,4),',',&
&meshgom(i,5),',',&
&meshgom(i,6),',',&
&meshgom(i,7),',',&
&meshgom(i,8)
end do
nyousosuugom=2*nx3*ny1*(nk+1)*nz
!載荷接点、拘束接点の設定
!片持ち
print'("*NSET,NSET=kousoku")'
do i=1,nitimen
nkousoku=i
print'(i6,a)',nkousoku,','
end do
print'("*NSET,NSET=saika")'
nsaikasettensuu=0
do i=1,nitimen
nsaika=nz*nitimen+i
print'(i6,a)',nsaika,','
nsaikasettensuu=nsaikasettensuu+1
end do
! 境界条件
print'("*BOUNDARY")'
print'("kousoku,1,3")'
! 木材の設定
print'("*MATERIAL,NAME=ELmoku")'
print'("*ELASTIC,TYPE=ORTHOROPIC")'
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
print'(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),','
print'(f12.4,a,a)',e66(4,4),',','293'
print'("*MATERIAL,NAME=ELkou")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*MATERIAL,NAME=ELgom")'
print'("*ELASTIC")'
print'("3.0d0, 0.45")'
print'("*SOLID SECTION,ELSET=Emoku,MATERIAL=ELmoku")'
print'("*SOLID SECTION,ELSET=Ekou,MATERIAL=ELkou")'
print'("*SOLID SECTION,ELSET=Egom,MATERIAL=ELgom")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe17.9)','saika',',',2,',',p/nsaikasettensuu
! datのオプション
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
print'("*NODE FILE")'
print'("U")'
print'("*END STEP")'
!
!
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



