!6節点要素のダイヤカット円筒を作るプログラム
implicit real*8 (a-h, o-z)
dimension x(10000),y(10000),z(10000),mesh6kaku(10000,7)&
&,mesh3kaku(10000,6),narabi(10000,10000)
pi=2.d0*asin(1.d0)!円周率
rentou=26.5d-3
atusa=1.00d-3
m=10!六角形の数(横方向)絶対偶数(4の倍数 2で割っても偶数の数)
n=45!六角形の数(z方向)絶対奇数
r=(rentou)/(2.0d0*m*sin(pi/(2.0d0*m)))!半径
b=(2.0d0*r*sin(pi/(2.0d0*m)))/(sqrt(3.0d0))!六角形の一辺の長さ
ell=((n-1)/2)*(2.0d0*b)+((n+1)/2)*b+b+b!軸方向長さ
m1=m*2
n1=n*2
ny=0 !三角形の分割の細かさ(4**nyだけ分割される)
th=pi/m!一辺が作る中心角
!奇数段の節点座標の入力(xとy)
!遇数段の節点座標の入力(xとy)
!節点座標の入力(z)
!節点番号と節点座標の入力と出力(少ない方の点下)
do j=1,n-1
do i=1,m
nset=i+(j-1)*((m+1)*2+m*2)
x(nset)=r*cos((i-1)*th+th/2.0d0)
y(nset)=r*sin((i-1)*th+th/2.0d0)
z(nset)=(j-1)*3.0d0*b+b/2.0d0
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!節点番号と節点座標の入力と出力(少ない方の点上)
do j=1,n-1
do i=1,m
nset=m+(m+1)*2+i+(j-1)*((m+1)*2+m*2)
x(nset)=r*cos((i-1)*th+th/2.0d0)
y(nset)=r*sin((i-1)*th+th/2.0d0)
z(nset)=(j-1)*3.0d0*b+2.0d0*b+b/2.0d0
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!節点番号と節点座標の入力と出力(多い点の方下)
do j=1,n-1
do i=1,m+1
nset=m+i+(j-1)*(2*m+2*(m+1))
x(nset)=r*cos((i-1)*th)
y(nset)=r*sin((i-1)*th)
z(nset)=b/2.0d0+(j-1)*3.0d0*b+b/2.0d0
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!節点番号と節点座標の入力と出力(多い点の方上)
do j=1,n-1
do i=1,m+1
nset=m+m+1+i+(j-1)*(2*m+2*(m+1))
x(nset)=r*cos((i-1)*th)
y(nset)=r*sin((i-1)*th)
z(nset)=b/2.0d0+b+(j-1)*3.0d0*b+b/2.0d0
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!六角形のグループの作成
!奇数段(六角形がm個の段)
!print'("*ELEMENT, TYPE=S6, ELSET=Eall")'
do j=1,n-1
do i=1,m
nyou=i+(j-1)*(2*m-1)
mesh6kaku(nyou,1)=i+(j-1)*(2*(m+1)+2*m)
mesh6kaku(nyou,2)=mesh6kaku(nyou,1)+m
mesh6kaku(nyou,3)=mesh6kaku(nyou,2)+1
mesh6kaku(nyou,4)=mesh6kaku(nyou,2)+m+1
mesh6kaku(nyou,5)=mesh6kaku(nyou,4)+1
mesh6kaku(nyou,6)=mesh6kaku(nyou,5)+m
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou&
!&,',',mesh6kaku(nyou,1)&
!&,',',mesh6kaku(nyou,2)&
!&,',',mesh6kaku(nyou,3)&
!&,',',mesh6kaku(nyou,4)&
!&,',',mesh6kaku(nyou,5)&
!&,',',mesh6kaku(nyou,6)
end do
end do
!奇数段(六角形がm個の段)
do j=1,n-2
do i=1,m-1
nyou=m+i+(j-1)*(2*m-1)
mesh6kaku(nyou,1)=(m+1)+m+1+i+(j-1)*(2*(m+1)+2*m)
mesh6kaku(nyou,2)=mesh6kaku(nyou,1)+m
mesh6kaku(nyou,3)=mesh6kaku(nyou,2)+1
mesh6kaku(nyou,4)=mesh6kaku(nyou,2)+m
mesh6kaku(nyou,5)=mesh6kaku(nyou,4)+1
mesh6kaku(nyou,6)=mesh6kaku(nyou,5)+m
!print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
!&nyou&
!&,',',mesh6kaku(nyou,1)&
!&,',',mesh6kaku(nyou,2)&
!&,',',mesh6kaku(nyou,3)&
!&,',',mesh6kaku(nyou,4)&
!&,',',mesh6kaku(nyou,5)&
!&,',',mesh6kaku(nyou,6)
end do
end do
!中心の節点をmesh6kakuにいれる。
n6kakusuu=m*(n+1)/2+(m-1)*(n-1)/2
nsettensuu=m*(n+1)+(m+1)*(n+1)
do i=1,n6kakusuu
x(i+nsettensuu)=(x(mesh6kaku(i,3))+x(mesh6kaku(i,4)))/2.0d0
y(i+nsettensuu)=(y(mesh6kaku(i,3))+y(mesh6kaku(i,4)))/2.0d0
z(i+nsettensuu)=(z(mesh6kaku(i,3))+z(mesh6kaku(i,4)))/2.0d0
mesh6kaku(i,7)=i+nsettensuu
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&i+nsettensuu&
!&,',',x(i+nsettensuu)&
!&,',',y(i+nsettensuu)&
!&,',',z(i+nsettensuu)
end do
!do i=1,n6kakusuu
!print'(i6,a,i6)',&
!&i,',',mesh6kaku(i,7)
!end do
!六角形のを6個の三角形に分ける
do i=1,n6kakusuu
mesh3kaku((i-1)*6+1,1)=mesh6kaku(i,1)
mesh3kaku((i-1)*6+1,2)=mesh6kaku(i,3)
mesh3kaku((i-1)*6+1,3)=mesh6kaku(i,7)
mesh3kaku((i-1)*6+2,1)=mesh6kaku(i,3)
mesh3kaku((i-1)*6+2,2)=mesh6kaku(i,5)
mesh3kaku((i-1)*6+2,3)=mesh6kaku(i,7)
mesh3kaku((i-1)*6+3,1)=mesh6kaku(i,5)
mesh3kaku((i-1)*6+3,2)=mesh6kaku(i,6)
mesh3kaku((i-1)*6+3,3)=mesh6kaku(i,7)
mesh3kaku((i-1)*6+4,1)=mesh6kaku(i,6)
mesh3kaku((i-1)*6+4,2)=mesh6kaku(i,4)
mesh3kaku((i-1)*6+4,3)=mesh6kaku(i,7)
mesh3kaku((i-1)*6+5,1)=mesh6kaku(i,4)
mesh3kaku((i-1)*6+5,2)=mesh6kaku(i,2)
mesh3kaku((i-1)*6+5,3)=mesh6kaku(i,7)
mesh3kaku((i-1)*6+6,1)=mesh6kaku(i,2)
mesh3kaku((i-1)*6+6,2)=mesh6kaku(i,1)
mesh3kaku((i-1)*6+6,3)=mesh6kaku(i,7)
end do
!do i=1,n6kakusuu*6
!print'(i6,a,i6,a,i6,a,i6)',&
!&i&
!&,',',mesh3kaku(i,1)&
!&,',',mesh3kaku(i,2)&
!&,',',mesh3kaku(i,3)
!end do
!中間節点を入れて6節点要素にする
do i=1,nsettensuu+n6kakusuu*7+((n+1)/2)*(m*2*2+m+1)&
&+((n-1)/2)*m+(2*m+1)*2+2*(n-1)/2
do j=1,nsettensuu+n6kakusuu*7+((n+1)/2)*(m*2*2+m+1)&
&+((n-1)/2)*m+(2*m+1)*2+2*(n-1)/2
narabi(i,j)=-1
end do
end do
iset=nsettensuu+n6kakusuu
do i=1,n6kakusuu*6
!1-2
if(narabi(mesh3kaku(i,1),mesh3kaku(i,2))==-1) then
iset=iset+1
narabi(mesh3kaku(i,1),mesh3kaku(i,2))=iset !頂点1-2間と
narabi(mesh3kaku(i,2),mesh3kaku(i,1))=iset !頂点2-1間とに節点番号isetをふり
mesh3kaku(i,4)=iset
x(iset)=(x(mesh3kaku(i,1))+x(mesh3kaku(i,2)))/2.d0
y(iset)=(y(mesh3kaku(i,1))+y(mesh3kaku(i,2)))/2.d0
z(iset)=(z(mesh3kaku(i,1))+z(mesh3kaku(i,2)))/2.d0
else
mesh3kaku(i,4)=narabi(mesh3kaku(i,1),mesh3kaku(i,2))
end if
!2-3
if(narabi(mesh3kaku(i,2),mesh3kaku(i,3))==-1) then
iset=iset+1
narabi(mesh3kaku(i,2),mesh3kaku(i,3))=iset !頂点2-3間と
narabi(mesh3kaku(i,3),mesh3kaku(i,2))=iset !頂点3-2間とに節点番号isetをふり
mesh3kaku(i,5)=iset
x(iset)=(x(mesh3kaku(i,2))+x(mesh3kaku(i,3)))/2.d0
y(iset)=(y(mesh3kaku(i,2))+y(mesh3kaku(i,3)))/2.d0
z(iset)=(z(mesh3kaku(i,2))+z(mesh3kaku(i,3)))/2.d0
else
mesh3kaku(i,5)=narabi(mesh3kaku(i,2),mesh3kaku(i,3))
end if
!3-1
if(narabi(mesh3kaku(i,3),mesh3kaku(i,1))==-1) then
iset=iset+1
narabi(mesh3kaku(i,3),mesh3kaku(i,1))=iset !頂点3-1間と
narabi(mesh3kaku(i,1),mesh3kaku(i,3))=iset !頂点1-3間とに節点番号isetをふり
mesh3kaku(i,6)=iset
x(iset)=(x(mesh3kaku(i,1))+x(mesh3kaku(i,3)))/2.d0
y(iset)=(y(mesh3kaku(i,1))+y(mesh3kaku(i,3)))/2.d0
z(iset)=(z(mesh3kaku(i,1))+z(mesh3kaku(i,3)))/2.d0
else
mesh3kaku(i,6)=narabi(mesh3kaku(i,1),mesh3kaku(i,3))
end if
end do
!端の節点を入れる
xmainasu=r*cos((-1)*th+th/2.0d0)
ymainasu=r*sin((-1)*th+th/2.0d0)
xhajihidari=(x(1)+xmainasu)/2.0d0
xhajimigi=(-1)*xhajihidari
yhajihidari=(y(1)+ymainasu)/2.0d0
yhajimigi=(-1)*yhajihidari
zhajihidari=0.00d0
zhajimigi=0.00d0
!端
do k=1,2
do i=1,2
nset=iset+(i-1)*(2*m)+1+(k-1)*(m*2+1+((n-1)/2)*2)
x(nset)=((-1)**(i+1))*(xhajihidari)
y(nset)=((-1)**(i+1))*(yhajihidari)
z(nset)=(k-1)*ell
print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
do k=1,(n-1)/2
do i=1,2
nset=iset+2*m+1+i+(k-1)*2
x(nset)=((-1)**(i+1))*(xhajihidari)
y(nset)=((-1)**(i+1))*(yhajihidari)
z(nset)=(k-1)*3.0d0*b+3.0d0*b
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
do k=1,2
do i=1,m-1
nset=iset+2+1+(i-1)*2+(k-1)*(m*2+1+((n-1)/2)*2)
x(nset)=x(nsettensuu+m+i)
y(nset)=y(nsettensuu+m+i)
z(nset)=(k-1)*ell
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
do k=1,2
do i=1,m
nset=iset+1+1+(i-1)*2+(k-1)*(m*2+1+((n-1)/2)*2)
x(nset)=x(i)
y(nset)=y(i)
z(nset)=(k-1)*ell
!print'(i6,a,f12.8,a,f12.8,a,f12.8)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
!端の要素を入れる
!左端
nyou=n6kakusuu*6+1
mesh3kaku(nyou,1)=iset+2
mesh3kaku(nyou,2)=mesh6kaku(1,1)
mesh3kaku(nyou,3)=iset+1
mesh3kaku(nyou+1,1)=mesh6kaku(1,1)
mesh3kaku(nyou+1,2)=mesh6kaku(1,2)
mesh3kaku(nyou+1,3)=iset+1
!do i=1,m*4*2+2*(n-1)/2
!print'(i6,a,i6,a,i6,a,i6)',&
!&i+n6kakusuu*6&
!&,',',mesh3kaku(i+n6kakusuu*6,1)&
!&,',',mesh3kaku(i+n6kakusuu*6,2)&
!&,',',mesh3kaku(i+n6kakusuu*6,3)
!end do
!右端
nyou=n6kakusuu*6+4*m-2+1
mesh3kaku(nyou,1)=mesh6kaku(m,3)
mesh3kaku(nyou,2)=mesh6kaku(m,1)
mesh3kaku(nyou,3)=iset+2*m+1
mesh3kaku(nyou+1,1)=mesh6kaku(m,1)
mesh3kaku(nyou+1,2)=iset+2*m
mesh3kaku(nyou+1,3)=iset+2*m+1
!左端先端
nyou=n6kakusuu*6+4*m+6*(n-1)/2+1
mesh3kaku(nyou,1)=mesh6kaku(n6kakusuu-m+1,4)
mesh3kaku(nyou,2)=mesh6kaku(n6kakusuu-m+1,6)
mesh3kaku(nyou,3)=iset+2*m+1+2*(n-1)/2+1
mesh3kaku(nyou+1,1)=mesh6kaku(n6kakusuu-m+1,6)
mesh3kaku(nyou+1,2)=iset+2*m+1+2*(n-1)/2+2
mesh3kaku(nyou+1,3)=iset+2*m+1+2*(n-1)/2+1
!右端先端
nyou=n6kakusuu*6+4*m+6*(n-1)/2+4*m-2+1
mesh3kaku(nyou,1)=iset+2*m+1+2*(n-1)/2+2*m
mesh3kaku(nyou,2)=mesh6kaku(n6kakusuu,6)
mesh3kaku(nyou,3)=iset+2*m+1+2*(n-1)/2+2*m+1
mesh3kaku(nyou+1,1)=mesh6kaku(n6kakusuu,6)
mesh3kaku(nyou+1,2)=mesh6kaku(n6kakusuu,5)
mesh3kaku(nyou+1,3)=iset+2*m+1+2*(n-1)/2+2*m+1
!
!do i=1,m*4*2+6*(n-1)/2
!print'(i6,a,i6,a,i6,a,i6)',&
!&i+n6kakusuu*6&
!&,',',mesh3kaku(i+n6kakusuu*6,1)&
!&,',',mesh3kaku(i+n6kakusuu*6,2)&
!&,',',mesh3kaku(i+n6kakusuu*6,3)
!end do
!
!端の中頃(z=0)
do i=1,m-1
mesh3kaku((i-1)*4+1+n6kakusuu*6+2,1)=iset+2+2+(i-1)*2
mesh3kaku((i-1)*4+1+n6kakusuu*6+2,2)=mesh6kaku(i+1,1)
mesh3kaku((i-1)*4+1+n6kakusuu*6+2,3)=iset+2+1+(i-1)*2
mesh3kaku((i-1)*4+2+n6kakusuu*6+2,1)=mesh6kaku(i+1,1)
mesh3kaku((i-1)*4+2+n6kakusuu*6+2,2)=mesh6kaku(i+1,2)
mesh3kaku((i-1)*4+2+n6kakusuu*6+2,3)=iset+2+1+(i-1)*2
mesh3kaku((i-1)*4+3+n6kakusuu*6+2,1)=mesh6kaku(i,3)
mesh3kaku((i-1)*4+3+n6kakusuu*6+2,2)=mesh6kaku(i,1)
mesh3kaku((i-1)*4+3+n6kakusuu*6+2,3)=iset+2+1+(i-1)*2
mesh3kaku((i-1)*4+4+n6kakusuu*6+2,1)=mesh6kaku(i,1)
mesh3kaku((i-1)*4+4+n6kakusuu*6+2,2)=iset+2+(i-1)*2
mesh3kaku((i-1)*4+4+n6kakusuu*6+2,3)=iset+2+1+(i-1)*2
end do
!先端の中頃(z=0)
do i=1,m-1
mesh3kaku((i-1)*4+1+n6kakusuu*6+4*m+6*(n-1)/2+2,1)=mesh6kaku(n6kakusuu-m+1+i,4)
mesh3kaku((i-1)*4+1+n6kakusuu*6+4*m+6*(n-1)/2+2,2)=mesh6kaku(n6kakusuu-m+1+i,6)
mesh3kaku((i-1)*4+1+n6kakusuu*6+4*m+6*(n-1)/2+2,3)=iset+2*m+1+2*(n-1)/2+2+1+(i-1)*2
mesh3kaku((i-1)*4+2+n6kakusuu*6+4*m+6*(n-1)/2+2,1)=mesh6kaku(n6kakusuu-m+1+i,6)
mesh3kaku((i-1)*4+2+n6kakusuu*6+4*m+6*(n-1)/2+2,2)=iset+2*m+1+2*(n-1)/2+2+2+(i-1)*2
mesh3kaku((i-1)*4+2+n6kakusuu*6+4*m+6*(n-1)/2+2,3)=iset+2*m+1+2*(n-1)/2+2+1+(i-1)*2
mesh3kaku((i-1)*4+3+n6kakusuu*6+4*m+6*(n-1)/2+2,1)=iset+2*m+1+2*(n-1)/2+2+(i-1)*2
mesh3kaku((i-1)*4+3+n6kakusuu*6+4*m+6*(n-1)/2+2,2)=mesh6kaku(n6kakusuu-m+i,6)
mesh3kaku((i-1)*4+3+n6kakusuu*6+4*m+6*(n-1)/2+2,3)=iset+2*m+1+2*(n-1)/2+2+1+(i-1)*2
mesh3kaku((i-1)*4+4+n6kakusuu*6+4*m+6*(n-1)/2+2,1)=mesh6kaku(n6kakusuu-m+i,6)
mesh3kaku((i-1)*4+4+n6kakusuu*6+4*m+6*(n-1)/2+2,2)=mesh6kaku(n6kakusuu-m+i,5)
mesh3kaku((i-1)*4+4+n6kakusuu*6+4*m+6*(n-1)/2+2,3)=iset+2*m+1+2*(n-1)/2+2+1+(i-1)*2
end do
!間の端左
do i=1,(n-1)/2
mesh3kaku(1+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(1+(2*m-1)*(i-1),4)
mesh3kaku(1+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(1+(2*m-1)*(i-1),6)
mesh3kaku(1+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+1+(i-1)*2
mesh3kaku(2+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(1+(2*m-1)*(i-1),6)
mesh3kaku(2+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(m+1+(2*m-1)*(i-1),4)
mesh3kaku(2+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+1+(i-1)*2
mesh3kaku(3+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(1+(2*m-1)*i,1)
mesh3kaku(3+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(1+(2*m-1)*i,2)
mesh3kaku(3+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+1+(i-1)*2
end do
!間の端右
do i=1,(n-1)/2
mesh3kaku(1+3+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(2*m+m-1+(2*m-1)*(i-1),3)
mesh3kaku(1+3+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(m+m-1+(2*m-1)*(i-1),5)
mesh3kaku(1+3+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+2+(i-1)*2
mesh3kaku(2+3+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(m+m-1+(2*m-1)*(i-1),5)
mesh3kaku(2+3+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(m+m-1+(2*m-1)*(i-1),3)
mesh3kaku(2+3+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+2+(i-1)*2
mesh3kaku(3+3+n6kakusuu*6+4*m+(i-1)*6,1)=mesh6kaku(m+(2*m-1)*(i-1),6)
mesh3kaku(3+3+n6kakusuu*6+4*m+(i-1)*6,2)=mesh6kaku(m+(2*m-1)*(i-1),5)
mesh3kaku(3+3+n6kakusuu*6+4*m+(i-1)*6,3)=iset+2*m+1+2+(i-1)*2
end do
!
!do i=1,m*4*2+6*(n-1)/2
!print'(i6,a,i6,a,i6,a,i6)',&
!&i+n6kakusuu*6&
!&,',',mesh3kaku(i+n6kakusuu*6,1)&
!&,',',mesh3kaku(i+n6kakusuu*6,2)&
!&,',',mesh3kaku(i+n6kakusuu*6,3)
!end do
!
!中間節点を入れて6節点要素にする
do i=1,(2*m+1)*2+2*(n-1)/2
do j=1,(2*m+1)*2+2*(n-1)/2
narabi(i+iset,j+iset)=-1
end do
end do
iset=iset+(2*m+1)*2+2*(n-1)/2
do i=1,4*m*2+6*(n-1)/2
!1-2
if(narabi(mesh3kaku(i+n6kakusuu*6,1),mesh3kaku(i+n6kakusuu*6,2))==-1) then
iset=iset+1
narabi(mesh3kaku(i+n6kakusuu*6,1),mesh3kaku(i+n6kakusuu*6,2))=iset !頂点1-2間と
narabi(mesh3kaku(i+n6kakusuu*6,2),mesh3kaku(i+n6kakusuu*6,1))=iset !頂点2-1間とに節点番号isetをふり
mesh3kaku(i+n6kakusuu*6,4)=iset
x(iset)=(x(mesh3kaku(i+n6kakusuu*6,1))+x(mesh3kaku(i+n6kakusuu*6,2)))/2.d0
y(iset)=(y(mesh3kaku(i+n6kakusuu*6,1))+y(mesh3kaku(i+n6kakusuu*6,2)))/2.d0
z(iset)=(z(mesh3kaku(i+n6kakusuu*6,1))+z(mesh3kaku(i+n6kakusuu*6,2)))/2.d0
else
mesh3kaku(i+n6kakusuu*6,4)=narabi(mesh3kaku(i+n6kakusuu*6,1),mesh3kaku(i+n6kakusuu*6,2))
end if
!2-3
if(narabi(mesh3kaku(i+n6kakusuu*6,2),mesh3kaku(i+n6kakusuu*6,3))==-1) then
iset=iset+1
narabi(mesh3kaku(i+n6kakusuu*6,2),mesh3kaku(i+n6kakusuu*6,3))=iset !頂点2-3間と
narabi(mesh3kaku(i+n6kakusuu*6,3),mesh3kaku(i+n6kakusuu*6,2))=iset !頂点3-2間とに節点番号isetをふり
mesh3kaku(i+n6kakusuu*6,5)=iset
x(iset)=(x(mesh3kaku(i+n6kakusuu*6,2))+x(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
y(iset)=(y(mesh3kaku(i+n6kakusuu*6,2))+y(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
z(iset)=(z(mesh3kaku(i+n6kakusuu*6,2))+z(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
else
mesh3kaku(i+n6kakusuu*6,5)=narabi(mesh3kaku(i+n6kakusuu*6,2),mesh3kaku(i+n6kakusuu*6,3))
end if
!3-1
if(narabi(mesh3kaku(i+n6kakusuu*6,3),mesh3kaku(i+n6kakusuu*6,1))==-1) then
iset=iset+1
narabi(mesh3kaku(i+n6kakusuu*6,3),mesh3kaku(i+n6kakusuu*6,1))=iset !頂点3-1間と
narabi(mesh3kaku(i+n6kakusuu*6,1),mesh3kaku(i+n6kakusuu*6,3))=iset !頂点1-3間とに節点番号isetをふり
mesh3kaku(i+n6kakusuu*6,6)=iset
x(iset)=(x(mesh3kaku(i+n6kakusuu*6,1))+x(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
y(iset)=(y(mesh3kaku(i+n6kakusuu*6,1))+y(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
z(iset)=(z(mesh3kaku(i+n6kakusuu*6,1))+z(mesh3kaku(i+n6kakusuu*6,3)))/2.d0
else
mesh3kaku(i+n6kakusuu*6,6)=narabi(mesh3kaku(i+n6kakusuu*6,1),mesh3kaku(i+n6kakusuu*6,3))
end if
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=S6, ELSET=Eall")'
do i=1,n6kakusuu*6+4*m*2+6*(n-1)/2
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i&
&,',',mesh3kaku(i,1)&
&,',',mesh3kaku(i,2)&
&,',',mesh3kaku(i,3)&
&,',',mesh3kaku(i,4)&
&,',',mesh3kaku(i,5)&
&,',',mesh3kaku(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,3")'
!材料の設定
print'("*MATERIAL,NAME=EL")'
print'("*ELASTIC")'
print'("69000.0,0.30000")'
!要素の設定,どの材料を使うか等
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,',',-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
!

