!6節点要素のダイヤカット円筒を作るプログラム
implicit real*8 (a-h, o-z)
dimension x(200000),y(200000),z(200000),mesh(200000,6)&
&,mesh2(200000,6),narabi(20000,20000)
character filename*128,jobnakami*88
ny=1 !三角形の分割の細かさ(4**nyだけ分割される)
open(8,file='jobdaiya1')
do m=4,20,2!周方向分割数を2ずつ変化させる
do n=4,20,2!高さ方向分割数2ずつ変化させる
write(jobnakami, '("ccx_1.7 d", i2.2, "-", i2.2, "-", i2.2,&
&">d", i2.2, "-", i2.2, "-", i2.2,".kekka")')m,n,2**ny&
&,m,n,2**ny
write(8,*) jobnakami
end do
end do
do m=4,20,2!周方向分割数を2ずつ変化させる
do n=4,20,2!高さ方向分割数2ずつ変化させる
write(filename, '("d", i2.2, "-", i2.2, "-", i2.2,".inp")')&
&m,n,2**ny 
open (7,file=filename, status='replace')
pi=2.d0*asin(1.d0)!円周率
b=53d-3*pi
h=104.5d-3 !高さ
atusa=0.8d-3
r=b/(2*m*sin(pi/m))-atusa/(2*cos(pi/m))!外半径。OFFSET=0を考慮。
m1=m*2
n1=n*2
nyouiti=2*m !一段の要素数
th=2.d0*pi/m!一辺が作る中心角
d=(1-cos(pi/m))*(b/(2*m*sin(pi/m)))
dh=sqrt((h/n)**2-d**2)!一段あたりの鉛直高さ
dh2=dh*n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)
!奇数段の節点座標の入力(xとy)
!遇数段の節点座標の入力(xとy)
!節点座標の入力(z)
!節点番号と節点座標の入力と出力(奇数段)
do j=0,n/2
do i=1,m
nset=i+j*2*m
x(nset)=r*cos((i-1)*th)
y(nset)=r*sin((i-1)*th)
z(nset)=dh*j*2
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nsetki,',',x(i),',',y(i),',',z(1+j*2)
end do
end do
!節点番号と節点座標の入力と出力(偶数段)
do j=0,n/2-1
do i=1,m
nset=i+j*2*m+m
x(nset)=r*cos((i-1)*th+th/2.0)
y(nset)=r*sin((i-1)*th+th/2.0)
z(nset)=dh+dh*j*2
!print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
!&nsetgu,',',x(i+m),',',y(i+m),',',z(2+j*2)
end do
end do
!中間節点の節点座標の入力と出力(凹み部奇数段)
do j=0,n/2
do i=1,m
nset=i+j*(2*m+2*m*2)+m*(n+1)
x(nset)=r*(cos((i-1)*th)+cos(i*th))/2.0
y(nset)=r*(sin((i-1)*th)+sin(i*th))/2.0
z(nset)=dh*j*2
!print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
!&nsetkiheko,',',x(i+m*2),',',y(i+m*2),',',z(2+j*2)
end do
end do
!中間節点の節点座標の入力と出力(凹み部遇数段)
do j=0,n/2-1
do i=1,m
nset=i+j*(2*m+2*m*2)+m*(n+1)+m+2*m
x(nset)=r*(cos((i-1)*th+th/2.0)+cos(i*th+th/2.0))/2.0
y(nset)=r*(sin((i-1)*th+th/2.0)+sin(i*th+th/2.0))/2.0
z(nset)=dh*j*2+dh
!print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
!&nsetguheko,',',x(i+m*3),',',y(i+m*3),',',z(2+j*2)
end do
end do
!中間節点の節点座標の入力と出力(段の中間節点)
r1=(r*(cos(0.)+cos(th/2.0))/2.0)/cos(th/4.0)!中間節点の半径
do j=0,n-1
do i=1,m*2
nset=i+j*(2*m+m)+m*(n+1)+m
x(nset)=r1*(cos((i-1)*th/2.0+th/4.0))
y(nset)=r1*(sin((i-1)*th/2.0+th/4.0))
z(nset)=dh*j+dh/2.0
!print'(i6,a,f9.6,a,f9.6,a,f9.6)',&
!&nsetguheko,',',x(i+m*3),',',y(i+m*3),',',z(2+j*2)
end do
end do
!三角形の要素を更に分割するための節点番号に-1を入れる
nset=m*2*(n*2+1)
do i=1,nset
do j=1,nset
narabi(i,j)=-1
end do
end do
!要素番号とその要素の3節点の節点番号の入力と出力
!逆三角形の要素(例外)奇数
!print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
do j=1,n/2
nyou1=nyouiti+(j-1)*2*nyouiti
mesh(nyou1,3)=1+m*(j-1)*2
mesh(nyou1,1)=mesh(nyou1,3)+m
mesh(nyou1,2)=mesh(nyou1,1)+m-1
mesh(nyou1,6)=m*(n+1)+m+1+(m+m*2)*(j-1)*2
mesh(nyou1,5)=m*(n+1)+m+m*2+(m+m*2)*(j-1)*2
mesh(nyou1,4)=m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou1,',',mesh(nyou1,1),',',mesh(nyou1,2),',',&
!&mesh(nyou1,3),',',mesh(nyou1,4),',',mesh(nyou1,5)&
!&,',',mesh(nyou1,6)
end do
!逆三角形の要素(例外)偶数
do j=1,n/2
nyou2=-1+j*2*nyouiti
mesh(nyou2,3)=m*2+m*(j-1)*2
mesh(nyou2,1)=mesh(nyou2,3)+1
mesh(nyou2,2)=mesh(nyou2,1)+m-1
mesh(nyou2,5)=m*(n+1)+m+m*2+m+m*2-1+(m+m*2)*(j-1)*2
mesh(nyou2,6)=mesh(nyou2,5)+1
mesh(nyou2,4)=mesh(nyou2,6)+m
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou2,',',mesh(nyou2,1),',',mesh(nyou2,2),',',&
!&mesh(nyou2,3),',',mesh(nyou2,4),',',mesh(nyou2,5)&
!&,',',mesh(nyou2,6)
end do
!三角形の要素(例外)奇数
do j=1,n/2
nyou3=nyouiti-1+(j-1)*2*nyouiti
mesh(nyou3,2)=1+m*(j-1)*2
mesh(nyou3,3)=mesh(nyou3,2)+m*2-1
mesh(nyou3,1)=mesh(nyou3,2)+m-1
mesh(nyou3,4)=m*(n+1)+m+(m+m*2)*(j-1)*2
mesh(nyou3,6)=m*(n+1)+m+m*2-1+(m+m*2)*(j-1)*2
mesh(nyou3,5)=mesh(nyou3,6)+1
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou3,',',mesh(nyou3,1),',',mesh(nyou3,2),',',&
!&mesh(nyou3,3),',',mesh(nyou3,4),',',mesh(nyou3,5)&
!&,',',mesh(nyou3,6)
end do
!三角形の要素(例外)偶数
do j=1,n/2
nyou4=j*2*nyouiti
mesh(nyou4,2)=m+1+m*(j-1)*2
mesh(nyou4,3)=mesh(nyou4,2)+m
mesh(nyou4,1)=mesh(nyou4,2)+m-1
mesh(nyou4,4)=m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2
mesh(nyou4,5)=m*(n+1)+m+m*2+m+1+(m+m*2)*(j-1)*2
mesh(nyou4,6)=m*(n+1)+m+m*2+m+m*2+(m+m*2)*(j-1)*2
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou4,',',mesh(nyou4,1),',',mesh(nyou4,2),',',&
!&mesh(nyou4,3),',',mesh(nyou4,4),',',mesh(nyou4,5)&
!&,',',mesh(nyou4,6)
end do
!三角形要素奇数
do j=1,n/2
do i=1,m-1
nyou5=-1+i*2+(j-1)*2*nyouiti
mesh(nyou5,1)=i+m*(j-1)*2
mesh(nyou5,2)=mesh(nyou5,1)+1
mesh(nyou5,3)=mesh(nyou5,1)+m
mesh(nyou5,4)=i+m*(n+1)+(m*2+m*2*2)*(j-1)
mesh(nyou5,6)=i*2-1+m*(n+1)+m+(m+m*2)*(j-1)*2
mesh(nyou5,5)=mesh(nyou5,6)+1
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',&
!&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)&
!&,',',mesh(nyou5,6)
end do
end do
!三角形要素遇数
do j=1,n/2
do i=1,m-1
nyou5=nyouiti+2+(i-1)*2+(j-1)*2*nyouiti
mesh(nyou5,1)=i+m*(j-1)*2+m
mesh(nyou5,2)=mesh(nyou5,1)+1
mesh(nyou5,3)=mesh(nyou5,1)+m+1
mesh(nyou5,4)=i+m*(n+1)+m+m*2+(m+m*2)*(j-1)*2
mesh(nyou5,6)=i*2+m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2
mesh(nyou5,5)=mesh(nyou5,6)+1
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',&
!&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)&
!&,',',mesh(nyou5,6)
end do
end do
!逆三角形要素奇数
do j=1,n/2
do i=1,m-1
nyou5=2+(i-1)*2+(j-1)*2*nyouiti
mesh(nyou5,3)=i+m*(j-1)*2+1
mesh(nyou5,1)=mesh(nyou5,3)+m
mesh(nyou5,2)=mesh(nyou5,3)+m-1
mesh(nyou5,5)=i*2+m*(n+1)+m+(m+m*2)*(j-1)*2
mesh(nyou5,6)=mesh(nyou5,5)+1
mesh(nyou5,4)=i+m*(n+1)+m+m*2+(m+m*2)*(j-1)*2
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',&
!&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)&
!&,',',mesh(nyou5,6)
end do
end do
!逆三角形要素遇数
do j=1,n/2
do i=1,m-1
nyou5=nyouiti+1+(i-1)*2+(j-1)*2*nyouiti
mesh(nyou5,3)=i+m*(j-1)*2+m
mesh(nyou5,1)=mesh(nyou5,3)+1+m
mesh(nyou5,2)=mesh(nyou5,3)+m
mesh(nyou5,5)=i*2-1+m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2
mesh(nyou5,6)=mesh(nyou5,5)+1
mesh(nyou5,4)=i+m*(n+1)+m+m*2+m+m*2+(m+m*2)*(j-1)*2
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',&
!&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)&
!&,',',mesh(nyou5,6)
end do
end do
!増やす時ためのmesh2(i,j)の設定
iset=nset
nyou=2*m*n
!mesh2にmeshを入れる
do i=1,nyou
do j=1,6
mesh2(i,j)=mesh(i,j)
end do 
end do
!
do 
if (ny==0) exit
!
do i=1,nyou
mesh2(1+(i-1)*4,1)=mesh(i,1)
mesh2(1+(i-1)*4,2)=mesh(i,4)
mesh2(1+(i-1)*4,3)=mesh(i,6)
mesh2(2+(i-1)*4,1)=mesh(i,4)
mesh2(2+(i-1)*4,2)=mesh(i,2)
mesh2(2+(i-1)*4,3)=mesh(i,5)
mesh2(3+(i-1)*4,1)=mesh(i,6)
mesh2(3+(i-1)*4,2)=mesh(i,5)
mesh2(3+(i-1)*4,3)=mesh(i,3)
mesh2(4+(i-1)*4,1)=mesh(i,5)
mesh2(4+(i-1)*4,2)=mesh(i,6)
mesh2(4+(i-1)*4,3)=mesh(i,4)
end do
!中間節点を入れて要素を増やす
do i=1,nyou
!1-4
if(narabi(mesh(i,1),mesh(i,4))==-1) then
iset=iset+1
narabi(mesh(i,1),mesh(i,4))=iset !頂点1-4間と
narabi(mesh(i,4),mesh(i,1))=iset !頂点4-1間とに節点番号isetをふり
mesh2(1+(i-1)*4,4)=iset
x(iset)=(x(mesh(i,1))+x(mesh(i,4)))/2.d0
y(iset)=(y(mesh(i,1))+y(mesh(i,4)))/2.d0
z(iset)=(z(mesh(i,1))+z(mesh(i,4)))/2.d0
else
mesh2(1+(i-1)*4,4)=narabi(mesh(i,1),mesh(i,4))
end if
!4-6
if(narabi(mesh(i,4),mesh(i,6))==-1) then
iset=iset+1
narabi(mesh(i,4),mesh(i,6))=iset !頂点1-4間と
narabi(mesh(i,6),mesh(i,4))=iset !頂点4-1間とに節点番号isetをふり
mesh2(1+(i-1)*4,5)=iset
mesh2(4+(i-1)*4,5)=iset
x(iset)=(x(mesh(i,4))+x(mesh(i,6)))/2.d0
y(iset)=(y(mesh(i,4))+y(mesh(i,6)))/2.d0
z(iset)=(z(mesh(i,4))+z(mesh(i,6)))/2.d0
else
mesh2(1+(i-1)*4,5)=narabi(mesh(i,4),mesh(i,6))
mesh2(4+(i-1)*4,5)=narabi(mesh(i,4),mesh(i,6))
end if
!1-6
if(narabi(mesh(i,1),mesh(i,6))==-1) then
iset=iset+1
narabi(mesh(i,1),mesh(i,6))=iset !頂点1-6間と
narabi(mesh(i,6),mesh(i,1))=iset !頂点6-1間とに節点番号isetをふり
mesh2(1+(i-1)*4,6)=iset
x(iset)=(x(mesh(i,1))+x(mesh(i,6)))/2.d0
y(iset)=(y(mesh(i,1))+y(mesh(i,6)))/2.d0
z(iset)=(z(mesh(i,1))+z(mesh(i,6)))/2.d0
else
mesh2(1+(i-1)*4,6)=narabi(mesh(i,1),mesh(i,6))
end if
!2-4
if(narabi(mesh(i,2),mesh(i,4))==-1) then
iset=iset+1
narabi(mesh(i,2),mesh(i,4))=iset !頂点1-4間と
narabi(mesh(i,4),mesh(i,2))=iset !頂点4-1間とに節点番号isetをふり
mesh2(2+(i-1)*4,4)=iset
x(iset)=(x(mesh(i,2))+x(mesh(i,4)))/2.d0
y(iset)=(y(mesh(i,2))+y(mesh(i,4)))/2.d0
z(iset)=(z(mesh(i,2))+z(mesh(i,4)))/2.d0
else
mesh2(2+(i-1)*4,4)=narabi(mesh(i,2),mesh(i,4))
end if
!2-5
if(narabi(mesh(i,2),mesh(i,5))==-1) then
iset=iset+1
narabi(mesh(i,2),mesh(i,5))=iset !頂点1-4間と
narabi(mesh(i,5),mesh(i,2))=iset !頂点4-1間とに節点番号isetをふり
mesh2(2+(i-1)*4,5)=iset
x(iset)=(x(mesh(i,2))+x(mesh(i,5)))/2.d0
y(iset)=(y(mesh(i,2))+y(mesh(i,5)))/2.d0
z(iset)=(z(mesh(i,2))+z(mesh(i,5)))/2.d0
else
mesh2(2+(i-1)*4,5)=narabi(mesh(i,2),mesh(i,5))
end if
!4-5
if(narabi(mesh(i,4),mesh(i,5))==-1) then
iset=iset+1
narabi(mesh(i,4),mesh(i,5))=iset !頂点1-4間と
narabi(mesh(i,5),mesh(i,4))=iset !頂点4-1間とに節点番号isetをふり
mesh2(2+(i-1)*4,6)=iset
mesh2(4+(i-1)*4,6)=iset
x(iset)=(x(mesh(i,4))+x(mesh(i,5)))/2.d0
y(iset)=(y(mesh(i,4))+y(mesh(i,5)))/2.d0
z(iset)=(z(mesh(i,4))+z(mesh(i,5)))/2.d0
else
mesh2(2+(i-1)*4,6)=narabi(mesh(i,4),mesh(i,5))
mesh2(4+(i-1)*4,6)=narabi(mesh(i,4),mesh(i,5))
end if
!5-6
if(narabi(mesh(i,5),mesh(i,6))==-1) then
iset=iset+1
narabi(mesh(i,5),mesh(i,6))=iset !頂点1-4間と
narabi(mesh(i,6),mesh(i,5))=iset !頂点4-1間とに節点番号isetをふり
mesh2(3+(i-1)*4,4)=iset
mesh2(4+(i-1)*4,4)=iset
x(iset)=(x(mesh(i,5))+x(mesh(i,6)))/2.d0
y(iset)=(y(mesh(i,5))+y(mesh(i,6)))/2.d0
z(iset)=(z(mesh(i,5))+z(mesh(i,6)))/2.d0
else
mesh2(3+(i-1)*4,4)=narabi(mesh(i,5),mesh(i,6))
mesh2(4+(i-1)*4,4)=narabi(mesh(i,5),mesh(i,6))
end if
!3-5
if(narabi(mesh(i,3),mesh(i,5))==-1) then
iset=iset+1
narabi(mesh(i,3),mesh(i,5))=iset !頂点1-4間と
narabi(mesh(i,5),mesh(i,3))=iset !頂点4-1間とに節点番号isetをふり
mesh2(3+(i-1)*4,5)=iset
x(iset)=(x(mesh(i,3))+x(mesh(i,5)))/2.d0
y(iset)=(y(mesh(i,3))+y(mesh(i,5)))/2.d0
z(iset)=(z(mesh(i,3))+z(mesh(i,5)))/2.d0
else
mesh2(3+(i-1)*4,5)=narabi(mesh(i,3),mesh(i,5))
end if
!3-6
if(narabi(mesh(i,3),mesh(i,6))==-1) then
iset=iset+1
narabi(mesh(i,3),mesh(i,6))=iset !頂点1-4間と
narabi(mesh(i,6),mesh(i,3))=iset !頂点4-1間とに節点番号isetをふり
mesh2(3+(i-1)*4,6)=iset
x(iset)=(x(mesh(i,3))+x(mesh(i,6)))/2.d0
y(iset)=(y(mesh(i,3))+y(mesh(i,6)))/2.d0
z(iset)=(z(mesh(i,3))+z(mesh(i,6)))/2.d0
else
mesh2(3+(i-1)*4,6)=narabi(mesh(i,3),mesh(i,6))
end if
end do
!mesh(i,j)をmesh2(i,j)に入れる
m1=m1*2
n1=n1*2
iset=m1*(n1+1)
do i=1,iset
do j=1,iset
narabi(i,j)=-1
end do
end do
nyou=nyou*4
do i=1,nyou
do j=1,6
mesh(i,j)=mesh2(i,j)
end do
end do
if (nyou==2*n*m*4**ny) exit
!
end do
!節点番号と座標
write(7,'("*NODE, NSET=Nall")')
do i=1,iset
write(7,'(i6,a,f12.8,a,f12.8,a,f12.8)')&
&i,',',x(i),',',y(i),',',z(i)
end do
!要素番号と節点番号
write(7,'("*ELEMENT, TYPE=S6, ELSET=Eall")')
do i=1,m*2*n*4**ny
write(7,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')&
&i,',',mesh2(i,1),',',mesh2(i,2),',',&
&mesh2(i,3),',',mesh2(i,4),',',mesh2(i,5)&
&,',',mesh2(i,6)
end do
!拘束節点と載荷節点の設定
write(7,'("*NSET, NSET=kotei")')
do i=1,iset
if (z(i)<0.00001.and.z(i)>-0.00001) then
write(7,'(i6,a)')i,','
end if
end do
write(7,'("*NSET, NSET=saika")')
do i=1,iset
if (z(i)<dh2+0.00001.and.z(i)>dh2-0.00001) then
write(7,'(i6,a)')i,','
end if
end do
!境界条件の設定
write(7,'("*BOUNDARY")')
write(7,'("kotei,1,3")')
write(7,'("saika,1,2")')
!材料の設定
write(7,'("*MATERIAL,NAME=EL")')
write(7,'("*ELASTIC")')
write(7,'("2843.4,0.313")')
!要素の設定,どの材料を使うか等
write(7,'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")')
write(7,'(1pe13.5)')atusa
!計算のオプション
write(7,'("*STEP")')
write(7,'("*STATIC,SOLVER=SPOOLES")')
!座屈解析のオプション
write(7,'("*BUCKLE")')
write(7,'("3,0.01")')
!荷重
write(7,'("*CLOAD")')
write(7,'(A,a,i6,a,1pe13.5)')'saika ',',',3,',',-2.00E-004/((m*2)*(2**ny))
!出力のオプション
write(7,'("*NODE PRINT,NSET=Nall")')
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")')
close(7)
end do
end do

end
!

