implicit real*8(a-h,o-z)
dimension x(20000),y(20000),z(20000),&
&mesh(20000,6),narabi(20000,20000)
!open(7,file='setten.inp')
!open(8,file='youso.inp')
!open(9,file='syoryou.inp')
a=43.3015d-3/2.0d0!格子の縦の半分のながさ
b=9.5238d-3/2.0d0!格子の横の半分のながさ
h=10.0d-3!たかさ
atusa=1.0d-3
p=5.0d-6
nz=4
nx=4
ny=10
ntatehan=9!縦板の数
nyokohan=2!横板の数
nyoko=ntatehan*2
ntate=nyokohan*2
pset=p/((ny+1)*ntatehan)
ell=a*ntate
xn=b/nx
zn=a/nz
yn=h/ny
!節点の入力
nxa=nyoko*nx
nza=ntate*nz
nitimen=(nxa+1)*nyokohan+(nza+1)*ntatehan-ntatehan*nyokohan
nsetten=nitimen*(ny+1)
!print'("*NODE, NSET=Nall")'
do j=1,ny+1
do l=1,2
do m=1,ntatehan
do k=1,nz
nset=k+(m-1)*nz+(l-1)*(nitimen-nz*ntatehan)+(j-1)*nitimen
x(nset)=b+(m-1)*b*2.0d0
y(nset)=(j-1)*yn
z(nset)=(k-1)*zn+(l-1)*(a*(ntate-1)+zn)
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!
!print'("*NODE, NSET=Nall")'
do j=1,ny+1
do l=1,nyokohan-1
do m=1,ntatehan
do k=1,nz*2-1
nset=k+(m-1)*(nz*2-1)+(j-1)*nitimen&
&+(l-1)*((nz*2-1)*ntatehan+nxa+1)&
&+(nz*ntatehan+nxa+1)
x(nset)=b+(m-1)*b*2.0d0
y(nset)=(j-1)*yn
z(nset)=k*zn+(l-1)*a*2.0d0+a
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!
!print'("*NODE, NSET=Nall")'
do j=1,ny+1
do l=1,nyokohan
do i=1,nxa+1
nset=i&
&+(l-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen+nz*ntatehan
x(nset)=(i-1)*xn
y(nset)=(j-1)*yn
z(nset)=(l-1)*a*2.0d0+a
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
!要素の設定
!横板三角要素
do j=1,ny
do m=1,nyokohan
do i=1,nxa
nyou=1+(i-1)*2+(m-1)*nxa*2*ny+(j-1)*nxa*2
mesh(nyou,1)=i+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
mesh(nyou,2)=i+1+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
mesh(nyou,3)=i+1+nitimen+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
end do
end do
end do
!横板三角要素
do j=1,ny
do m=1,nyokohan
do i=1,nxa
nyou=2+(i-1)*2+(m-1)*nxa*2*ny+(j-1)*nxa*2
mesh(nyou,3)=i+nz*ntatehan&
&+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
mesh(nyou,2)=i+nitimen+nz*ntatehan&
&+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
mesh(nyou,1)=i+1+nitimen+nz*ntatehan&
&+(m-1)*(nxa+1+(nz*2-1)*ntatehan)&
&+(j-1)*nitimen
end do
end do
end do
nyou=nxa*ny*nyokohan*2
!縦板はじ
do j=1,ny
do l=1,2
do m=1,ntatehan
do k=1,nz-1
nyou=nyou+1
mesh(nyou,2)=k+(m-1)*nz+(l-1)*(nitimen-nz*ntatehan)&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=mesh(nyou,2)+1
mesh(nyou+1,2)=mesh(nyou+1,1)+1
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
end do
end do
!縦板交差1-1
do j=1,ny
do m=1,ntatehan
nyou=nyou+1
mesh(nyou,2)=nz+(m-1)*nz&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=nz*ntatehan+nx+1+(m-1)*nx*2+(j-1)*nitimen
mesh(nyou+1,2)=mesh(nyou,1)+nitimen
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
!縦板交差1-2kara
do j=1,ny
do l=1,nyokohan-1
do m=1,ntatehan
nyou=nyou+1
mesh(nyou,2)=nz*ntatehan+nxa+1&
&+nz*2-1+(l-1)*((nz*2-1)*ntatehan+nxa+1)&
&+(m-1)*(nz*2-1)&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=nz*ntatehan+nxa+1+(nz*2-1)*ntatehan&
&+(m-1)*nx*2+(l-1)*((nz*2-1)*ntatehan+nxa+1)&
&+nx+1+(j-1)*nitimen
mesh(nyou+1,2)=mesh(nyou,1)+nitimen
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
end do
!縦板交差2
do j=1,ny
do l=1,nyokohan-1
do m=1,ntatehan
nyou=nyou+1
mesh(nyou,2)=nz*ntatehan+nx+1+(m-1)*nx*2&
&+(l-1)*((nz*2-1)*ntatehan+nxa+1)&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=nz*ntatehan+nxa+1+1+(m-1)*(nz*2-1)&
&+(l-1)*((nz*2-1)*ntatehan+nxa+1)+(j-1)*nitimen
mesh(nyou+1,2)=mesh(nyou,1)+nitimen
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
end do
!縦板交差2-last
do j=1,ny
do m=1,ntatehan
nyou=nyou+1
mesh(nyou,2)=nz*ntatehan+nx+1&
&+((nz*2-1)*ntatehan+nxa+1)*(nyokohan-1)+(m-1)*nx*2&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=nz*ntatehan+nxa+1+1&
&+((nz*2-1)*ntatehan+nxa+1)*(nyokohan-1)&
&+(m-1)*nz&
&+(j-1)*nitimen
mesh(nyou+1,2)=mesh(nyou,1)+nitimen
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
!縦板間
do j=1,ny
do l=1,nyokohan-1
do m=1,ntatehan
do k=1,nz*2-2
nyou=nyou+1
mesh(nyou,2)=k+nz*ntatehan+nxa+1&
&+(m-1)*(nz*2-1)+(l-1)*((nz*2-1)*ntatehan+nxa+1)&
&+(j-1)*nitimen
mesh(nyou+1,1)=mesh(nyou,2)+nitimen
mesh(nyou,1)=k+1+nz*ntatehan+nxa+1&
&+(m-1)*(nz*2-1)+(l-1)*((nz*2-1)*ntatehan+nxa+1)+(j-1)*nitimen
mesh(nyou+1,2)=mesh(nyou,1)+nitimen
mesh(nyou,3)=mesh(nyou,2)+nitimen
mesh(nyou+1,3)=mesh(nyou,1)
nyou=nyou+1
end do
end do
end do
end do
!write(7,*) nsetten
!do i=1,nsetten
!write(7,'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)')&
!&i,',',x(i),',',y(i),',',z(i)
!end do
!write(8,*) nyou
!do i=1,nyou
!write(8,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')&
!&i,',',&
!&mesh(i,1),',',&
!&mesh(i,2),',',&
!&mesh(i,3)
!end do
!write(9,*) atusa
!write(9,*) p/(2*(ny+1)*ntatehan)
!write(9,*) e
!close(7)
!close(8)
!close(9)
!節点を増やしてs6要素にする
iset=nsetten
!narabiに-1をいれる
do j=1,iset
do i=1,iset
narabi(i,j)=-1
end do
end do
!中間節点を入れて要素を増やす
do i=1,nyou
!1-2
if(narabi(mesh(i,1),mesh(i,2))==-1) then
iset=iset+1
narabi(mesh(i,1),mesh(i,2))=iset !頂点1-2間と
narabi(mesh(i,2),mesh(i,1))=iset !頂点2-1間とに節点番号isetをふり
mesh(i,4)=iset
x(iset)=(x(mesh(i,1))+x(mesh(i,2)))/2.d0
y(iset)=(y(mesh(i,1))+y(mesh(i,2)))/2.d0
z(iset)=(z(mesh(i,1))+z(mesh(i,2)))/2.d0
else
mesh(i,4)=narabi(mesh(i,1),mesh(i,2))
end if
!2-3
if(narabi(mesh(i,2),mesh(i,3))==-1) then
iset=iset+1
narabi(mesh(i,2),mesh(i,3))=iset !頂点2-3間と
narabi(mesh(i,3),mesh(i,2))=iset !頂点3-2間とに節点番号isetをふり
mesh(i,5)=iset
x(iset)=(x(mesh(i,2))+x(mesh(i,3)))/2.d0
y(iset)=(y(mesh(i,2))+y(mesh(i,3)))/2.d0
z(iset)=(z(mesh(i,2))+z(mesh(i,3)))/2.d0
else
mesh(i,5)=narabi(mesh(i,2),mesh(i,3))
end if
!1-3
if(narabi(mesh(i,1),mesh(i,3))==-1) then
iset=iset+1
narabi(mesh(i,1),mesh(i,3))=iset !頂点1-3間と
narabi(mesh(i,3),mesh(i,1))=iset !頂点3-1間とに節点番号isetをふり
mesh(i,6)=iset
x(iset)=(x(mesh(i,1))+x(mesh(i,3)))/2.d0
y(iset)=(y(mesh(i,1))+y(mesh(i,3)))/2.d0
z(iset)=(z(mesh(i,1))+z(mesh(i,3)))/2.d0
else
mesh(i,6)=narabi(mesh(i,1),mesh(i,3))
end if
end do
!
print'("*NODE, NSET=Nall")'
do i=1,iset
print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
&i,',',x(i),',',y(i),',',z(i)
end do
print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
do i=1,nyou
print'(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)
end do
print'("*NSET, NSET=kotei")'
do i=1,iset
if (z(i)<0.00001.and.z(i)>-0.00001) then
print '(i6,a)',i,','
end if
end do
print'("*NSET, NSET=saika")'
do i=1,iset
if (z(i)<ell+0.00001.and.z(i)>ell-0.00001) then
print '(i6,a)',i,','
end if
end do
!境界条件の設定
print'("*BOUNDARY")'
print'("kotei,1,6")'
!材料の設定
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC")'
print'("3011.7,0.3290")'
!要素の設定,どの材料を使うか等
print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")'
print'(1pe13.5)',atusa
!計算のオプション
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
!荷重
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika ',',',2,',',p
!出力のオプション
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
print'("*EL PRINT,ELSET=Eall")'
print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE,POSITION=AVERAGED AT NODES")'
print'("S")'
print'("*END STEP")'

!
end 
