!一枚の板からダイヤカット円筒の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 !外側の高さ変更 fai=dh/(h/n)!角度 r3=r2+atusa/fai!高さを合わせるための半径 !print*,r !print*,r2 !print*,r3 !奇数段の節点座標の入力(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 !内側w 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !2014/9/5坪井が作った部分★ミ !ゼロ弾目上書き do i=1,m nset=i !外側 x(nset)=r3*cos((i-1)*th) y(nset)=r3*sin((i-1)*th) z(nset)=dage end do !最後の段上書き do j=n/2,n/2 do i=1,m nset=i+j*2*m !外側 x(nset)=r3*cos((i-1)*th) y(nset)=r3*sin((i-1)*th) z(nset)=dage+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