!一枚の板からダイヤカット円筒のstlファイルを作るプログラム
!角度による厚さの変化を考慮したもの
!内側の円筒を上に上げる
implicit real*8 (a-h, o-z)
dimension x(2000),y(2000),z(2000)&
&,houx(2000),houy(2000),houz(2000),houookisa(2500),&
&mesh(2000,3),mesh2(1000,3),mesh3(4000,4),&
&meshface(2500,3)
pi=2.d0*asin(1.d0)!円周率
b=53d-3*pi     !元となる長方形の横幅
h=104.5d-3    !元となる長方形の高さ
m=12 !周方向分割数
n=10 !高さ方向分割数
r=b/(2*m*sin(pi/m))!外半径
d=(1-cos(pi/m))*r
nyouiti=2*m !一段の要素数
th=2.d0*pi/m!一辺が作る中心角
dh=sqrt((h/n)**2-d**2)!一段あたりの鉛直高さ
atusa=0.8d-3!厚さ
r2=r-(atusa/cos(pi/m))!内側の半径
dh2=dh*r2/r
!内側の円筒を上に上げるdage
dage=((dh-dh2)*n)/2
!!
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)=dage+dh2*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)=dage+dh2+dh2*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
!face要素の設定(外側と内側)
do i=1,nyousosuu
do j=1,3
meshface(i,j)=mesh(i,j)
meshface(i+nyousosuu,j)=mesh2(i,j)
end do 
end do
!!face要素の設定(xy平面の上と下)
do i=1,m*2-1,2
meshface(nyousosuu*2+1+(i-1),1)=mesh(i,1)
meshface(nyousosuu*2+1+(i-1),2)=mesh(i,2)
meshface(nyousosuu*2+1+(i-1),3)=mesh2(i,1)
meshface(nyousosuu*2+2+(i-1),1)=mesh2(i,2)
meshface(nyousosuu*2+2+(i-1),2)=mesh2(i,1)
meshface(nyousosuu*2+2+(i-1),3)=mesh(i,2)
meshface(nyousosuu*2+1+m*2+(i-1),1)=mesh(i+nyousosuu-m*2,1)
meshface(nyousosuu*2+1+m*2+(i-1),2)=mesh(i+nyousosuu-m*2,2)
meshface(nyousosuu*2+1+m*2+(i-1),3)=mesh2(i+nyousosuu-m*2,1)
meshface(nyousosuu*2+2+m*2+(i-1),1)=mesh2(i+nyousosuu-m*2,2)
meshface(nyousosuu*2+2+m*2+(i-1),2)=mesh2(i+nyousosuu-m*2,1)
meshface(nyousosuu*2+2+m*2+(i-1),3)=mesh(i+nyousosuu-m*2,2)
end do
do i=1,nyousosuu
houx(i)=(y(mesh(i,2))-y(mesh(i,1)))*(z(mesh(i,3))-z(mesh(i,1)))&
&-(z(mesh(i,2))-z(mesh(i,1)))*(y(mesh(i,3))-y(mesh(i,1)))
houy(i)=(z(mesh(i,2))-z(mesh(i,1)))*(x(mesh(i,3))-x(mesh(i,1)))&
&-(x(mesh(i,2))-x(mesh(i,1)))*(z(mesh(i,3))-z(mesh(i,1)))
houz(i)=(x(mesh(i,2))-x(mesh(i,1)))*(y(mesh(i,3))-y(mesh(i,1)))&
&-(y(mesh(i,2))-y(mesh(i,1)))*(x(mesh(i,3))-x(mesh(i,1)))
end do
do i=1,nyousosuu
houx(i+nyousosuu)=houx(i)*(-1.0d0)
houy(i+nyousosuu)=houy(i)*(-1.0d0)
houz(i+nyousosuu)=houz(i)*(-1.0d0)
end do
do i=1+nyousosuu*2,nyousosuu*2+m*2
houx(i)=0.000d0
houy(i)=-1.000d0
houz(i)=0.0000d0
end do
do i=1+nyousosuu*2,nyousosuu*2+m*2
houx(i+m*2)=0.000d0
houy(i+m*2)=1.000d0
houz(i+m*2)=0.000d0
end do
do i=1,nyousosuu*2+m*2*2
houookisa(i)=sqrt(houx(i)**2+houy(i)**2+houz(i)**2)
end do
print '("solid")'
do i=1,nyousosuu
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,' facet normal ',houx(i)/houookisa(i),' ',houy(i)/houookisa(i)&
&,' ',houz(i)/houookisa(i)
print '("   outer loop")'
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,1)),' ',y(meshface(i,1)),' ',z(meshface(i,1))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,2)),' ',y(meshface(i,2)),' ',z(meshface(i,2))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,3)),' ',y(meshface(i,3)),' ',z(meshface(i,3))
print '("   endloop")'
print '(" endfacet")'
end do
do i=1+nyousosuu,nyousosuu*2
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,' facet normal ',houx(i)/houookisa(i),' ',houy(i)/houookisa(i)&
&,' ',houz(i)/houookisa(i)
print '("   outer loop")'
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,1)),' ',y(meshface(i,1)),' ',z(meshface(i,1))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,3)),' ',y(meshface(i,3)),' ',z(meshface(i,3))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,2)),' ',y(meshface(i,2)),' ',z(meshface(i,2))
print '("   endloop")'
print '(" endfacet")'
end do
do i=1+nyousosuu*2,nyousosuu*2+m*2
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,' facet normal ',houx(i)/houookisa(i),' ',houy(i)/houookisa(i)&
&,' ',houz(i)/houookisa(i)
print '("   outer loop")'
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,1)),' ',y(meshface(i,1)),' ',z(meshface(i,1))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,3)),' ',y(meshface(i,3)),' ',z(meshface(i,3))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,2)),' ',y(meshface(i,2)),' ',z(meshface(i,2))
print '("   endloop")'
print '(" endfacet")'
end do
do i=1+nyousosuu*2+m*2,nyousosuu*2+m*2*2
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,' facet normal ',houx(i)/houookisa(i),' ',houy(i)/houookisa(i)&
&,' ',houz(i)/houookisa(i)
print '("   outer loop")'
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,1)),' ',y(meshface(i,1)),' ',z(meshface(i,1))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,3)),' ',y(meshface(i,3)),' ',z(meshface(i,3))
print '(a,1pe13.5,a,1pe13.5,a,1pe13.5)'&
&,'     vertex  ',x(meshface(i,2)),' ',y(meshface(i,2)),' ',z(meshface(i,2))
print '("   endloop")'
print '(" endfacet")'
end do
print '("endsolid")'
end 
