!mentat用の3節点要素のダイヤカット円筒を作るプログラム
implicit real*8 (a-h, o-z)
dimension x(100000),y(100000),z(100000),mesh(1000000,3)&
&,mesh2(1000000,3),narabi(10000,10000)
pi=2.d0*asin(1.d0)!円周率
r=0.03 !半径
h=0.5d0 !高さ
m=4 !周方向分割数
n=4 !高さ方向分割数
m1=m
n1=n
ny=2
nyouiti=2*m !一段の要素数
th=2.d0*pi/m!一辺が作る中心角
dh=h/n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)
atusa=2.000d-4
!奇数段の節点座標の入力(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
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
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
!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
!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
!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
!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
!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
!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
!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
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3)
end do
end do
!mesh2にmeshを入れる
nyou=2*m*n
do i=1,nyou
do j=1,3
mesh2(i,j)=mesh(i,j)
end do
end do
!三角形の要素を更に分割するための節点番号に-1を入れる
nset=m*(n+1)
do i=1,nset
do j=1,nset
narabi(i,j)=-1
end do
end do
!増やす時のためのmesh2(i,j)の設定
iset=nset
nyou=2*m*n
do 
if (ny==0) exit
do i=1,nyou
mesh2(1+(i-1)*4,1)=mesh(i,1)
mesh2(2+(i-1)*4,2)=mesh(i,2)
mesh2(3+(i-1)*4,3)=mesh(i,3)
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をふり
mesh2(1+(i-1)*4,2)=iset
mesh2(2+(i-1)*4,1)=iset
mesh2(4+(i-1)*4,3)=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
mesh2(1+(i-1)*4,2)=narabi(mesh(i,1),mesh(i,2))
mesh2(2+(i-1)*4,1)=narabi(mesh(i,1),mesh(i,2))
mesh2(4+(i-1)*4,3)=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をふり
mesh2(2+(i-1)*4,3)=iset
mesh2(3+(i-1)*4,2)=iset
mesh2(4+(i-1)*4,1)=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
mesh2(2+(i-1)*4,3)=narabi(mesh(i,2),mesh(i,3))
mesh2(3+(i-1)*4,2)=narabi(mesh(i,2),mesh(i,3))
mesh2(4+(i-1)*4,1)=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をふり
mesh2(1+(i-1)*4,3)=iset
mesh2(3+(i-1)*4,1)=iset
mesh2(4+(i-1)*4,2)=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
mesh2(1+(i-1)*4,3)=narabi(mesh(i,1),mesh(i,3))
mesh2(3+(i-1)*4,1)=narabi(mesh(i,1),mesh(i,3))
mesh2(4+(i-1)*4,2)=narabi(mesh(i,1),mesh(i,3))
end if
end do
!isetに現在の節点数を入れる
!narabiに-1を入れなおす
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
!mesh2(i,j)をmesh2(i,j)に入れる
nyou=nyou*4
do i=1,nyou
do j=1,3
mesh(i,j)=mesh2(i,j)
end do
end do
if (nyou==m*2*n*4**ny) exit
end do
!節点の出力
print'("*NODE, NSET=Nall")'
do i=1,iset
print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
&i,',',x(i),',',y(i),',',z(i)
end do
!要素の出力
print'("*ELEMENT, TYPE=S3R, ELSET=Eall")'
do i=1,nyou
print'(i6,a,i6,a,i6,a,i6)',&
&i,',',mesh2(i,1),',',mesh2(i,2),',',mesh2(i,3)
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)>dh-0.00001.and.z(i)<dh+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'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")'
!shellの厚さ
print'(1pe13.5)',atusa
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika ',',',2,',',1.00E-003/m
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
!

