!等圧縮を受ける板
implicit real*8(a-h,o-z)
dimension x(4000000),y(4000000),z(4000000),e66(6,6),b(3,3)&
&,mesh(3000000,8)
character filename*128
!プログラムはメートル、MPa、
b1=9.0d-3!厚さ
h=5.d-1!高さ(m) (a)
ell=5d-1!幅(m)(b)
!p=1.d0
p=0.5d0
!p=2.d0
!p=1.d0!変更
!xyz方向の要素分割数
nx=10
ny=60
nz=60
!1つあたりの要素の長さ?
xn=b1/nx
yn=h/ny
zn=ell/nz
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
nitimen=(nx+1)*(ny+1)
do k=0,nz!z
do j=0,ny!
do i=0,nx!x方向分割数を担当
nset=1+i+j*(nx+1)+k*((nx+1)*(ny+1))
x(nset)=i*xn
y(nset)=j*yn
z(nset)=k*zn
end do
end do
end do 
!i5~6が整数,aが文字,1pe13.5が小数.1peとかはEの＊乗
!接点番号と座標の出力
nsettenzenbu=(nx+1)*(ny+1)*(nz+1)
print'("*NODE, NSET=Nall")'
do i=1,nsettenzenbu
print'(i8,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
    &i,',',x(i),',',y(i),',',z(i)
nsettensuu=i
end do
!要素番号とその要素の8接点の節点番号
nyou=0
do k=0,nz-1
do j=0,ny-1
do i=1,nx
nyou=nyou+1
mesh(nyou,1)=i+j*(nx+1)+k*((nx+1)*(ny+1))
mesh(nyou,2)=mesh(nyou,1)+(nx+1)
mesh(nyou,3)=mesh(nyou,2)+(nx+1)*(ny+1)
mesh(nyou,4)=mesh(nyou,1)+(nx+1)*(ny+1)
mesh(nyou,5)=mesh(nyou,1)+1
mesh(nyou,6)=mesh(nyou,5)+nx+1
mesh(nyou,7)=mesh(nyou,6)+(nx+1)*(ny+1)
mesh(nyou,8)=mesh(nyou,5)+(nx+1)*(ny+1)
end do
end do
end do
!要素の出力
print'("*ELEMENT, TYPE=C3D8, ELSET=E1")'
do i=1,nyou
print'(i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8,a,i8)',&
    &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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!載荷接点、拘束接点の設定
print'("*NSET,NSET=kousoku")'
do m=0,1
do l=0,1
do j=1,(ny/2)
do i=1,nx+1
nkousoku=i+(j-1)*(nx+1)+l*(nx+1)*((ny/2)+1)+m*nitimen*nz
print'(i8,a)',nkousoku
end do
end do
end do
end do
print'("*NSET,NSET=kousoku2")'
do l=0,1
do k=1,nz/2
do m=0,1
do i=1,nx+1
nkousoku2=i+m*(nx+1)*ny+(k-1)*nitimen+l*nitimen*((nz/2)+1)
print'(i8,a)',nkousoku2
end do
end do
end do
end do
print'("*NSET,NSET=kousoku3")'
do m=0,1
do i=1,nx+1
nkousoku3=(nx+1)*(ny/2)+i+m*nitimen*nz
print'(i8,a)',nkousoku3
end do
end do
print'("*NSET,NSET=kousoku4")'
!do m=0,1
do i=1,nx+1
!nkousoku4=nitimen*(nz/2)+i+m*(nx+1)*ny
nkousoku4=nitimen*(nz/2)+i
print'(i8,a)',nkousoku4
!end do
end do
!print'("*NSET,NSET=kousoku5")'
!do i=1,nsettensuu
!if ( h-0.00001d0<y(i).and. h+0.00001d0>y(i)) then
!nkousoku5=i
!print'(i8,a)',nkousoku5
!end if
!end do
!print'("*NSET,NSET=kousoku3")'
!do m=0,1
!nkousoku3=(nx+1)*(ny/2)+(nx/2)+1+m*nitimen*nz
!print'(i8,a)',nkousoku3
!end do
!print'("*NSET,NSET=kousoku4")'
!do m=0,1
!nkousoku4=nitimen*(nz/2)+(nx/2)+1+m*(nx+1)*ny
!print'(i8,a)',nkousoku4
!end do
print'("*NSET,NSET=kousokumen")'
do k=0,nz-2
do i=1,nx+1
nkmen=nitimen+(nx+1)*(ny/2)+i+k*nitimen
print'(i8,a)',nkmen
end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
print'("*NSET,NSET=saika")'
nsaiten=0
do k=0,nz
do i=1,nx+1
nsaika=i+k*nitimen
print'(i8,a)',nsaika
nsaiten=nsaiten+1
end do
end do
print'("*NSET,NSET=saika2")'
nsaiten2=0
do k=0,nz
do i=1,nx+1
nsaika2=nitimen-(nx+1)+i+k*nitimen
print'(i8,a)',nsaika2
nsaiten2=nsaiten2+1
end do
end do
!境界条件
print'("*BOUNDARY")'
print'("kousoku,3")'
print'("kousoku3,1,3")'
print'("kousokumen,2")'
print'("kousoku,1")'
!print'("kousoku4,1,3")'
!print'("kousoku2,1,2")'
! 木材の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*BUCKLE")'
print'("5,0.01")'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
print'("*CLOAD")'
print'(a6,a,i4,a,1pe13.5)','saika',',',2,',',p/nsaiten
print'(a6,a,i4,a,1pe13.5)','saika2',',',2,',',-p/nsaiten2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! datのオプション
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
print'("*EL PRINT,ELSET=E1")'
print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'
end
