!mentat用の3節点要素のダイヤカット円筒を作るプログラム
implicit real*8 (a-h, o-z)
dimension x(100000),y(100000),z(100000)&
&,x2(100000),y2(100000),z2(100000)&
&,mesh(30000,3),mesh2(30000,6),narabi(10000,10000)
pi=2.d0*asin(1.d0)!円周率
r=0.03 !半径
h=0.5d0 !高さ
m=100 !周方向分割数
n=20 !高さ方向分割数
th=2.d0*pi/m!一辺が作る中心角mentat
th2=2.d0*pi/(m*2)!一辺が作る中心角ccx
dh=h/n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)mentat
dh2=h/(n*2) !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)ccx
!奇数段の節点座標の入力(xとy)
!遇数段の節点座標の入力(xとy)
!節点座標の入力(z)
!節点番号と節点座標の入力と出力
open(7,file='kakutoumentat.inp')
open(8,file='kakutouccx.inp')
!write(7,'("*NODE, NSET=Nall")')
nset=0
do j=0,n
do i=1,m
nset=nset+1
x(nset)=r*cos((i-1)*th)
y(nset)=r*sin((i-1)*th)
z(nset)=dh*j*2
!write(7,'(i6,a,f12.8,a,f12.8,a,f12.8)')&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!write(8,'("*NODE, NSET=Nall")')
!write(8,'(i6,a,f12.8,a,f12.8,a,f12.8)')&
!&nset,',',x2(nset),',',y2(nset),',',z2(nset)
!要素番号とその要素の3節点の節点番号の設定(mentat)
!逆三角形の要素(例外)
!print'("*ELEMENT, TYPE=S3R, ELSET=Eall")'
do j=1,n
nyou=m*2*j
mesh(nyou,3)=1+m*(j-1)
mesh(nyou,1)=1+m*j
mesh(nyou,2)=m*(j+1)
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou,',',mesh(nyou,1),',',mesh(nyou,2),',',mesh(nyou,3)
end do
!三角形の要素(例外)
do j=1,n
nyou=m*2*j-1
mesh(nyou,2)=1+m*(j-1)
mesh(nyou,1)=m*j
mesh(nyou,3)=m*(j+1)
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou,',',mesh(nyou,1),',',mesh(nyou,2),',',mesh(nyou,3)
end do
!三角形要素
do j=1,n
do i=1,m-1
nyou=-1+i*2+(j-1)*2*m
mesh(nyou,1)=i+m*(j-1)
mesh(nyou,2)=i+1+m*(j-1)
mesh(nyou,3)=i+m*j
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou,',',mesh(nyou,1),',',mesh(nyou,2),',',mesh(nyou,3)
end do
end do
!逆三角形要素
do j=1,n
do i=1,m-1
nyou=i*2+(j-1)*2*m
mesh(nyou,3)=i+1+m*(j-1)
mesh(nyou,1)=i+1+m*j
mesh(nyou,2)=i+m*j
!print'(i6,a,i6,a,i6,a,i6)',&
!&nyou,',',mesh(nyou,1),',',mesh(nyou,2),',',mesh(nyou,3)
end do
end do
!
!要素番号とその要素の3節点の節点番号の設定(ccx)
!mesh2(i,j)にmesh(i,j) j=1,3をいれる。
do i=1,2*m*n
do j=1,3
mesh2(i,j)=mesh(i,j)
end do
end do
!narabi(i,j)に-1をいれる
do j=1,m*(n+1)
do i=1,m*(n+1)
narabi(i,j)=-1
end do 
end do
!中間節点をふっていく
iset=m*(n+1)
do i=1,2*m*n
!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(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
mesh2(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をふり
mesh2(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
mesh2(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をふり
mesh2(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
mesh2(i,6)=narabi(mesh(i,1),mesh(i,3))
end if
end do
!mentat用ファイルの作成
!節点出力
write(7,'("*NODE, NSET=Nall")')
do i=1,m*(n+1)
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=S3R, ELSET=Eall")')
do i=1,m*2*n
write(7,'(i6,a,i6,a,i6,a,i6)')&
&i,',',mesh(i,1),',',mesh(i,2),',',mesh(i,3)
end do
!拘束節点出力
write(7,'("*NSET, NSET=kotei")')
do i=1,m
write(7,'(i6,a)')i,','
end do
!載荷節点出力
write(7,'("*NSET, NSET=saika")')
do i=1,m
write(7,'(i6,a)')i+m*n,','
end do
!境界条件出力
write(7,'("*BOUNDARY")')
write(7,'("kotei,1,3")')
!材料定数出力
write(7,'("*MATERIAL,NAME=EL")')
write(7,'("*ELASTIC")')
write(7,'("69000.0,0.30000")')
!要素の材料の設定
write(7,'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")')
write(7,'("2.0E-004")')
write(7,'("*STEP")')
write(7,'("*STATIC,SOLVER=SPOOLES")')
!載荷条件
write(7,'("*CLOAD")')
write(7,'(A,a,i6,a,1pe13.5)')'saika ',',',2,',',1.00E-003/m
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")')
!ccx用ファイルの作成
!節点出力
write(8,'("*NODE, NSET=Nall")')
do i=1,2*m*(2*n+1)
write(8,'(i6,a,f12.8,a,f12.8,a,f12.8)')&
&i,',',x(i),',',y(i),',',z(i)
end do
!要素出力
write(8,'("*ELEMENT, TYPE=S6, ELSET=Eall")')
do i=1,m*2*n
write(8,'(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(8,'("*NSET, NSET=kotei")')
do i=1,m*2
write(8,'(i6,a)')i,','
end do
!載荷節点出力
write(8,'("*NSET, NSET=saika")')
do i=1,m
write(8,'(i6,a)')-1+i*2+4*m*n,','
end do
!境界条件出力
write(8,'("*BOUNDARY")')
write(8,'("kotei,1,3")')
!材料定数出力
write(8,'("*MATERIAL,NAME=EL")')
write(8,'("*ELASTIC")')
write(8,'("69000.0,0.30000")')
!要素の材料の設定
write(8,'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")')
write(8,'("2.0E-004")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
!載荷条件
write(8,'("*CLOAD")')
write(8,'(A,a,i6,a,1pe13.5)')'saika ',',',2,',',1.00E-003/m
write(8,'("*NODE PRINT,NSET=Nall")')
write(8,'("U")')
write(8,'("*EL PRINT,ELSET=Eall")')
write(8,'("S")')
write(8,'("*NODE FILE")')
write(8,'("U")')
write(8,'("*EL FILE,POSITION=AVERAGED AT NODES")')
write(8,'("S")')
write(8,'("*END STEP")')


end
!

