!mentat用の3節点要素のダイヤカット円筒を作るプログラム
implicit real*8 (a-h, o-z)
dimension x(2000),y(2000),z(2000),&
&mesh(2000,3),mesh2(2000,3),mesh3(4000,4)
pi=2.d0*asin(1.d0)!円周率
r=0.03 !半径
h=0.5d0 !高さ
m=4 !周方向分割数
n=2 !高さ方向分割数
nyouiti=2*m !一段の要素数
th=2.d0*pi/m!一辺が作る中心角
dh=h/n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)
atusa=3.000d-3
r2=r-atusa
nsettensuu=m*(n+1)
nyousosuu=2*m*n
!奇数段の節点座標の入力(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
x(nset+nsettensuu)=r2*cos((i-1)*th)
y(nset+nsettensuu)=r2*sin((i-1)*th)
z(nset+nsettensuu)=dh*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)
y(nset)=r*sin((i-1)*th+th/2)
z(nset)=dh+dh*j*2
x(nset+nsettensuu)=r2*cos((i-1)*th+th/2)
y(nset+nsettensuu)=r2*sin((i-1)*th+th/2)
z(nset+nsettensuu)=dh+dh*j*2
end do
end do
!nset=0
!do j=1,n+1
!do i=1,m
!nset=nset+1
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
!end do
!end do
!要素番号とその要素の3節点の節点番号の入力と出力
!逆三角形の要素(例外)奇数
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
mesh2(nyou1,3)=1+m*(j-1)*2+nsettensuu
mesh2(nyou1,1)=mesh(nyou1,3)+m+nsettensuu
mesh2(nyou1,2)=mesh(nyou1,1)+m-1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou1,',',mesh(nyou1,1),',',mesh(nyou1,2),',',mesh(nyou1,3)
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
mesh2(nyou2,3)=m*2+m*(j-1)*2+nsettensuu
mesh2(nyou2,1)=mesh(nyou2,3)+1+nsettensuu
mesh2(nyou2,2)=mesh(nyou2,1)+m-1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou2,',',mesh(nyou2,1),',',mesh(nyou2,2),',',mesh(nyou2,3)
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
mesh2(nyou3,2)=1+m*(j-1)*2+nsettensuu
mesh2(nyou3,3)=mesh(nyou3,2)+m*2-1+nsettensuu
mesh2(nyou3,1)=mesh(nyou3,2)+m-1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou3,',',mesh(nyou3,1),',',mesh(nyou3,2),',',mesh(nyou3,3)
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
mesh2(nyou4,2)=m+1+m*(j-1)*2+nsettensuu
mesh2(nyou4,3)=mesh(nyou4,2)+m+nsettensuu
mesh2(nyou4,1)=mesh(nyou4,2)+m-1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou4,',',mesh(nyou4,1),',',mesh(nyou4,2),',',mesh(nyou4,3)
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
mesh2(nyou5,1)=i+m*(j-1)*2+nsettensuu
mesh2(nyou5,2)=mesh(nyou5,1)+1+nsettensuu
mesh2(nyou5,3)=mesh(nyou5,1)+m+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3)
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
mesh2(nyou5,1)=i+m*(j-1)*2+m+nsettensuu
mesh2(nyou5,2)=mesh(nyou5,1)+1+nsettensuu
mesh2(nyou5,3)=mesh(nyou5,1)+m+1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3)
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
mesh2(nyou5,3)=i+m*(j-1)*2+1+nsettensuu
mesh2(nyou5,1)=mesh(nyou5,3)+m+nsettensuu
mesh2(nyou5,2)=mesh(nyou5,3)+m-1+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3)
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
mesh2(nyou5,3)=i+m*(j-1)*2+m+nsettensuu
mesh2(nyou5,1)=mesh(nyou5,3)+1+m+nsettensuu
mesh2(nyou5,2)=mesh(nyou5,3)+m+nsettensuu
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3)
end do
end do
!c3d4要素の設定
!一つ目の要素
do i=1,nyousosuu
mesh3((i-1)*3+1,1)=mesh(i,3)
mesh3((i-1)*3+1,2)=mesh(i,2)
mesh3((i-1)*3+1,3)=mesh(i,1)
mesh3((i-1)*3+1,4)=mesh2(i,1)
end do
!２つ目の要素
do i=1,nyousosuu
mesh3((i-1)*3+2,1)=mesh(i,3)
mesh3((i-1)*3+2,2)=mesh(i,2)
mesh3((i-1)*3+2,3)=mesh2(i,1)
mesh3((i-1)*3+2,4)=mesh2(i,2)
end do
!3つ目の要素
do i=1,nyousosuu
mesh3((i-1)*3+3,1)=mesh(i,3)
mesh3((i-1)*3+3,2)=mesh2(i,2)
mesh3((i-1)*3+3,3)=mesh2(i,1)
mesh3((i-1)*3+3,4)=mesh2(i,3)
end do
!!inpファイルを書く
print'("*NODE, NSET=Nall")'
do i=1,nsettensuu*2
print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',i,',',x(i),',',y(i),',',z(i)
end do
print'("*ELEMENT,TYPE=C3D4, ELSET=Eall")'
do i=1,m*2*n*3
print'(i6,a,i6,a,i6,a,i6,a,i6)',i,',',mesh3(i,1),',',&
&mesh3(i,2),',',mesh3(i,3),',',mesh3(i,4)
end do
print'("*NSET, NSET=kotei")'
do i=1,nsettensuu*2
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,nsettensuu*2
if (z(i)<h+0.00001.and.z(i)>h-0.00001) then
print '(i6,a)',i,','
end if
end do
!境界条件の設定
print'("*BOUNDARY")'
print'("kotei,1,3")'
!材料の設定
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC")'
print'("69000.0,0.30000")'
!要素の設定,どの材料を使うか等
print'("SOLID SECTION,ELSET=Eall,MATERIAL=EL")'
!計算のオプション
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
!荷重
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika ',',',3,',',-1.00E-003/((m*2)*(2**ny))
!出力のオプション
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
