implicit real*8(a-h,o-z)
dimension x(300000),y(300000),z(300000),e66(6,6),b(3,3)&
&,mesh(3000000,8)
character filename*128,filename2*48
!プログラムはメートル、MPa、
ell=2.0d0
h2=-2.0d-3!高さ(m)
do i3=100,490,10
write(filename, '("hako", i4.4, ".inp")') i3!
open  (7,file=filename, status='replace')
!寸法
b1=1.0d-3!幅(m)
b2=8.0d-3!幅(m)
h1=1.0d-3!高さ(m)
h2=h2+1.0d-2!長さ
p=1.0d-3
!木材のとき
!ヤング率
!ezz=7.694d3
!exx=ezz/25.0d0!x軸方向のヤング率
!eyy=ezz/25.0d0!y軸方向のヤング率
!ezz2=206.0d3!z軸方向のヤング率
!ポアソン比の入力
!poixy=0.016d0
!poixz=poixy
!poiyx=poixy
!poiyz=poixy
!poizx=0.4d0
!せん断弾性係数の入力
!gxy=ezz/15.0d0
!gxz=gxy
!gyz=gxy
!xyz方向の要素分割数
nx1=2
nx2=6
nx3=2
ny1=2
ny2=6
ny3=2
nz=1000
!1つあたりの要素の長さ?
xn1=b1/nx1
xn2=b2/nx2
xn3=b1/nx3
yn1=h1/ny1
yn2=h2/ny2
yn3=h1/ny3
zn=ell/nz
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える
!左下
nitimen=((nx1+nx2+nx3+1)*(ny1+ny3+2)+(nx1+nx3+2)*(ny2-1))
nitidan1=(nx1+nx2+nx3+1)

do k=0,nz!z
do j=0,ny1!
do i=0,nx1!x方向分割数を担当
nset=1+i+j*nitidan1+k*nitimen
x(nset)=i*xn1            
y(nset)=j*yn1            
z(nset)=k*zn             
end do                   
end do                 
end do

!真ん中下
do k=0,nz!z
do j=0,ny1!
do i=0,nx2-2!x方向分割数を担当
nset=nx1+2+i+j*nitidan1+k*nitimen
x(nset)=b1+(1+i)*xn2
y(nset)=j*yn1
z(nset)=k*zn
end do
end do
end do

!右下                                                 
do k=0,nz!z
do j=0,ny1!
do i=0,nx3!x方向分割数を担当
nset=nx1+nx2+1+i+j*nitidan1+k*nitimen
x(nset)=b1+b2+i*xn3            
y(nset)=j*yn1            
z(nset)=k*zn             
end do                   
end do                 
end do

!左
nitidan2=nx1+nx3+2

do k=0,nz!z
do j=0,ny2-2!
do i=0,nx1!x方向分割数を担当
nset=((nx1+nx2+nx3+1)*(ny1+1))+1+i+j*nitidan2+k*nitimen
x(nset)=i*xn1
y(nset)=h1+((1+j)*yn2)
z(nset)=k*zn
end do
end do
end do

!右
do k=0,nz!z
do j=0,ny2-2!
do i=0,nx3!x方向分割数を担当
nset=((nx1+nx2+nx3+1)*(ny1+1))+nx1+2+i+j*nitidan2+k*nitimen
x(nset)=b1+b2+i*xn3
y(nset)=h1+((1+j)*yn2)
z(nset)=k*zn
end do
end do
end do

!左上                                                            
do k=0,nz!z
do j=0,ny3!
do i=0,nx1!x方向分割数を担当
nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+1+i+j*nitidan1+k*nitimen
x(nset)=i*xn1            
y(nset)=h1+h2+j*yn3            
z(nset)=k*zn             
end do                   
end do                 
end do                  

!真ん中上          
do k=0,nz!z
do j=0,ny3!
do i=0,nx2-2!x方向分割数を担当
nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+2+i+j*nitidan1+k*nitimen
x(nset)=b1+(1+i)*xn2
y(nset)=h1+h2+j*yn3
z(nset)=k*zn
end do
end do
end do

!右上
do k=0,nz!z
do j=0,ny3!
do i=0,nx3!x方向分割数を担当
nset=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+nx2+1+i+j*nitidan1+k*nitimen
x(nset)=b1+b2+i*xn3
y(nset)=h1+h2+j*yn3
z(nset)=k*zn
end do
end do
end do
                                                                                
!i5~6が整数,aが文字,f9,6が小数.1peとかはEの＊乗
!接点番号と座標の出力
nsettenzenbu=nitimen*(nz+1)
write(7,'("*NODE, NSET=Nall")')
do i=1,nsettenzenbu
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 k=0,nz-1
do j=0,ny1
do i=1,nx1
nyou=nyou+1
mesh(nyou,1)=i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!真ん中下
do k=0,nz-1
do j=0,ny1-1
do i=1,nx2
nyou=nyou+1
mesh(nyou,1)=nx1+i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!右下
do k=0,nz-1
do j=0,ny1-1
do i=1,nx3
nyou=nyou+1
mesh(nyou,1)=nx1+nx2+i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!左
do k=0,nz-1
do j=0,ny2-2
do i=1,nx1
nyou=nyou+1
mesh(nyou,1)=nitidan1*(ny1+1)+i+j*nitidan2+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan2
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan2
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!右
do k=0,nz-1
do j=0,ny2-2
do i=1,nx3
nyou=nyou+1
mesh(nyou,1)=nitidan1*ny1+nx1+nx2+i+j*nitidan2+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan2
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan2
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!左上
do k=0,nz-1
do j=0,ny3-1
do i=1,nx1
nyou=nyou+1
mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-1)+i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!真ん中上                                                
do k=0,nz-1
do j=0,ny3-1
do i=1,nx2
nyou=nyou+1
mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-1)+nx1+i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!右上
do k=0,nz-1
do j=0,ny3
do i=1,nx3
nyou=nyou+1
mesh(nyou,1)=nitidan1*(ny1+1)+nitidan2*(ny2-2)+nx1+1+i+j*nitidan1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+nitidan1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nitidan1
mesh(nyou,7)=mesh(nyou,6)+nitimen
mesh(nyou,8)=mesh(nyou,5)+nitimen
end do
end do
end do

!要素の出力
write(7,'("*ELEMENT, TYPE=C3D8, ELSET=E1")')
do i=1,nyou
write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)'),&
    &i,',',&
    &mesh(i,1),',',&
    &mesh(i,2),',',&
    &mesh(i,3),',',&
    &mesh(i,4),',',&
    &mesh(i,5),',',&
    &mesh(i,6),',',&
    &mesh(i,7),',',&
    &mesh(i,8)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!固定端の節点群を書きだす
write(7,'("*NSET,NSET=kotei")')
do i=1,nitimen
nkotei=i
write(7,'(i6,a)'),nkotei
end do
!載荷点の節点群を書きだす
write(7,'("*NSET,NSET=saika")')
do i=1,nitimen
nsaika=nitimen*nz+i
write(7,'(i6,a)'),nsaika
end do
!境界条件
write(7,'("*BOUNDARY")')
write(7,'("kotei,1,3")')
! 木材の設定
write(7,'("*MATERIAL,NAME=EL1")')
write(7,'("*ELASTIC")')
write(7,'("206.d3, 0.3")')
write(7,'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")')
! 載荷条件
write(7,'("*STEP")')
write(7,'("*STATIC,SOLVER=SPOOLES")')
!print'("*BUCKLE")'
!print'("5,0.01")'
write(7,'("*CLOAD")')
write(7,'(A,a,i6,a,1pe13.5)'),'saika',',',2,',',p/nitimen
! datのオプション
write(7,'("*NODE PRINT,NSET=saika")')
write(7,'("U")')
write(7,'("*EL PRINT,ELSET=Eall")')
write(7,'("S")')
write(7,'("*NODE FILE")')
write(7,'("U")')
write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")')
write(7,'("S")')
write(7,'("*END STEP")')
!
end do
!jobの設定
do i=100,490,10
write(filename2, '("ccx_2.5 hako", i4.4,">hako", i4.4,".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
!
!****** せん断補正係数を入力 ******
!write(7,*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!write(7,*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!write(7,*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!
!
return
end
