!onsi21220.f90改良版
!ただの箱断面
implicit real*8(a-h,o-z)
dimension x(1000000),y(1315500),z(1315500),e66(6,6),b(3,3)&
&,meshmoku(300000,8),meshkou(300000,8)
character filename*128,filename2*48
ell=2.5d0
b1=50.0d-2
b2=0.0
h1=0.0
h2=50.0d-2
hall=2*h1+h2
p=1.000d-3!荷重
!ここを変えたら下の!jobの設定 のi3も同じく変えて
do i3=100,210,10!1から100まで1刻み
write(filename, '("hako15-ell", i3.3, ".inp")') i3!
open  (7,file=filename, status='replace')
b2=b2+2.0d-2
b1=b1-4.0d-2
h1=h1+2.0d-2
h2=h2-4.0d-2
!xyz方向の要素分割数
nx1=10!木材
nx2=8!鋼板
ny1=10!上部分と下部分の高さ方向の分割数
ny2=8!間の鋼板の高さ方向の分割数(3以上)
nz=200
!ヤング率の入力
ezz=206.0d3
exx=ezz!x軸方向のヤング率
eyy=ezz!y軸方向のヤング率
ezz2=206.0d3!z軸方向のヤング率
!ポアソン比の入力
poixy=0.3d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.3d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/(2.0d0*(1.0d0+poixy))
gxz=gxy
gyz=gxy
!1つあたりの要素の長さ?
xn1=b1/nx1
xn2=b2/nx2
yn1=h1/ny1
yn2=h2/ny2
zn=ell/nz
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
!下と上部分の節点　鋼板!
do m=0,1
do k=0,nz
do j=0,ny1
do l=0,1
do i=0,nx2
nset=1+i+l*(nx2+nx1)+j*(nx2*2+nx1+1)&
&+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))&
&+m*((nx2*2+nx1+1)*(ny1+1)+(nx2+1)*2*(ny2-1))
!上と下あるのでmを動かして表す
!kは面の一面すべての節点を掛けてz方向に動かす

x(nset)=i*xn2+l*(b1+b2) !iが0~xn2までいってi=0,l=1になる
y(nset)=j*yn1+m*(h1+h2)
z(nset)=k*zn
end do
end do
end do 
end do
end do
!下と上部分の節点 木材
do m=0,1
do k=0,nz
do j=0,ny1
do i=0,nx1-2
nset=nx2+1+1+i+j*(nx2*2+nx1+1)&
&+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))&
&+m*((nx2*2+nx1+1)*(ny1+1)+(nx2+1)*2*(ny2-1))
!上と下2つあるのでmを動かして表す
x(nset)=i*xn1+b2+xn1 !iが0~xn2までいってi=0,l=1になる
y(nset)=j*yn1+m*(h1+h2)
z(nset)=k*zn
end do
end do
end do
end do

!真ん中部分の鋼材の節点
do k=0,nz
do j=0,ny2-2
do l=0,1
do i=0,nx2
nset=(1+nx1+2*nx2)*(ny1+1)+1+i+l*(1+nx2)+j*2*(nx2+1)&
&+k*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))
x(nset)=i*xn2+l*(b1+b2)
y(nset)=j*yn2+h1+yn2
z(nset)=k*zn
end do 
end do
end do
end do

!i5~6が整数,aが文字,f9,6が小数.1peとかはEの＊乗
!接点番号と座標の出力
write(7,'("*NODE, NSET=Nall")')
do i=1,(nz+1)*(2*(nx1+nx2*2+1)*(ny1+1)+2*(nx2+1)*(ny2-1))
write(7,'(i6,a,f9.6,a,f9.6,a,f9.6)'),&
&i,',',x(i),',',y(i),',',z(i)
nsettensuu=i
end do

!要素番号とその要素の8接点の接点番号
!下と上部分木材の要素
nyou=0
do m=0,1
do k=0,nz-1
do j=0,ny1-1
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=nx2+i+(2*nx2+nx1+1)*j&
&+k*((2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*((2*nx2+nx1+1)*(ny1+1)+2*(nx2+1)*(ny2-1))
meshmoku(nyou,2)=meshmoku(nyou,1)+(2*nx2+nx1+1)
meshmoku(nyou,3)=meshmoku(nyou,2)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshmoku(nyou,4)=meshmoku(nyou,1)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,5)+(2*nx2+nx1+1)
meshmoku(nyou,7)=meshmoku(nyou,6)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshmoku(nyou,8)=meshmoku(nyou,5)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
end do
end do
end do
end do
!上部分と下部分の鋼板
nyou=0
do l=0,1
do m=0,1
do k=0,nz-1 
do j=0,ny1-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=i+(2*nx2+nx1+1)*j&
&+k*((2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*((2*nx2+nx1+1)*(ny1+1)+2*(nx2+1)*(ny2-1))&
&+l*(nx2+nx1)
meshkou(nyou,2)=meshkou(nyou,1)+(2*nx2+nx1+1)
meshkou(nyou,3)=meshkou(nyou,2)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshkou(nyou,4)=meshkou(nyou,1)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,5)+(2*nx2+nx1+1)
meshkou(nyou,7)=meshkou(nyou,6)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshkou(nyou,8)=meshkou(nyou,5)+(2*nx2+nx1+1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
end do
end do 
end do
end do
end do
!真ん中の鋼板(最初の一段目)
do m=0,1
do k=0,nz-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=i+(nx2*2+nx1+1)*ny1&
&+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*(nx2+nx1)
meshkou(nyou,2)=(2*nx2+nx1+1)*(ny1+1)+i&
&+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*(nx2+1)
meshkou(nyou,3)=meshkou(nyou,2)+(nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
meshkou(nyou,4)=meshkou(nyou,1)+(nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1)
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 m=0,1
do k=0,nz-1
do j=0,ny2-3!ny2>3
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=(2*nx2+nx1+1)*(ny1+1)+i&
&+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*(nx2+1)+j*(nx2+1+nx2+1)
meshkou(nyou,2)=meshkou(nyou,1)+nx2+1+nx2+1
meshkou(nyou,3)=meshkou(nyou,2)+(nx2+1)*(ny2-1)*2+(nx2*2+1+nx1)*(ny1+1)*2
meshkou(nyou,4)=meshkou(nyou,1)+(nx2+1)*(ny2-1)*2+(nx2*2+1+nx1)*(ny1+1)*2
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 m=0,1
do k=0,nz-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=(nx2*2+nx1+1)*(ny1+1)+(nx2+1)*(ny2-2)*2+i&
&+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*(nx2+1)
meshkou(nyou,2)=(nx2*2+nx1+1)*(ny1+1)+(nx2+1)*(ny2-1)*2+i&
&+k*((nx2*2+1+nx1)*(ny1+1)*2+2*(nx2+1)*(ny2-1))&
&+m*(nx2+nx1)
meshkou(nyou,3)=meshkou(nyou,2)+(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2
meshkou(nyou,4)=meshkou(nyou,1)+(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2
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
write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")')
do i=1,2*nx2*(ny1*2+ny2)*nz
write(7,'(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
write(7,'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")')
do i=1,2*nx1*ny1*nz
write(7,'(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

!載荷接点、拘束接点の設定
write(7,'("*NSET,NSET=kousoku")')
do i=1,(nx2*2+nx1+1)*2*(ny1+1)+(nx2+1)*(ny2-1)*2
nkousoku=i
write(7,'(i6,a)'),nkousoku,','
end do
write(7,'("*NSET,NSET=saika")')
nsaikasettensuu=0
do i=1,nsettensuu    !節点数
if (ell-0.00002<z(i).and.z(i)<ell+0.00002)then
write(7,'(i6,a)'),i,','
nsaikasettensuu=nsaikasettensuu+1
end if
end do
write(7,'("*NSET,NSET=saikasita")')
do i=1,nsettensuu    !節点数
if (ell-0.00002<z(i).and.z(i)<ell+0.00002&
    &.and.y(i)>-0.00001d0.and.y(i)<0.00001d0&
        &.and.b2+0.0001d0<x(i).and.x(i)<b1+b2-0.0001d0)then
write(7,'(i6,a)'),i,','
end if
end do
! 境界条件
write(7,'("*BOUNDARY")')
write(7,'("kousoku,1,3")')
! 木材の設定
write(7,'("*MATERIAL,NAME=ELmoku")')
write(7,'("*ELASTIC,TYPE=ORTHOROPIC")')
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
write(7,'(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(7,'(f12.4,a,a)'),e66(4,4),',','293'
write(7,'("*MATERIAL,NAME=ELkou")')
write(7,'("*ELASTIC")')
write(7,'("206.d3, 0.3")')
write(7,'("*SOLID SECTION,ELSET=Emoku,MATERIAL=ELmoku")')
write(7,'("*SOLID SECTION,ELSET=Ekou,MATERIAL=ELkou")')
! 載荷条件
write(7,'("*STEP")')
write(7,'("*STATIC,SOLVER=SPOOLES")')
write(7,'("*CLOAD")')
write(7,'(A,a,i6,a,1pe13.5)'),'saika',',',2,',',p/nsaikasettensuu
! datのオプション
write(7,'("*NODE print,NSET=saika")')
write(7,'("U")')
write(7,'("*NODE FILE")')
write(7,'("U")')
write(7,'("*END STEP")')
!
close (7)
!
end do
!

!jobの設定
do i=100,210,10
write(filename2, '("ccx_2.5 hako15-ell", i3.3,">hako15-ell", i3.3,".kekka" )') i,i
open  (8,file='job')
write(8,*) filename2
end do
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
