! utf-8 ! 文字コードはutf-8 ! 自由形式のfortranで書いています。 ! g77でコンパイルする場合は、 !g77 -ffree-form -o smv smv.f ! のようにコンパイルしてください。 ! gfortran でコンパイルする場合は文字コードをutf-8にしてください。 implicit real*8 (a-h, o-z) dimension x(3000000),y(3000000),z(3000000),& & xu(3000000),yu(3000000),zu(3000000),& & mesh(1000000,6),narabi(10000,10000),& & mesh2(1000000,6) character mozi*20,mozi0*20 !!!!!!!!!!!! !!!!!!!!!!!! pi=2.d0*asin(1.d0) r=1.d0/2.d0/pi !半径 h=50.d-2 !高さ m=10 !周方向分割数 n=5 !高さ方向分割数 pi=2.d0*asin(1.d0) atusa=9.d-3 call daiya(r,h,m,n) r=r-atusa call daiya(r,h,m,n) close(8) !!!!!!!!!!!! !!!!!!!!!!!! !!!!!!!!!!!! open(9,file='daiya3d.obj') !yousosuu.fで数えた要素数ファイル read(9,*) mozi0 !1行目の文字を読み込む read(9,*) nset !節点数の読み込み do i=0,nset-1 read(9,*) x(i),y(i),z(i) !各節点の座標値の読み込み end do ! open(7,file='yousosuu.txt') !yousosuu.fで数えた要素数ファイル read(7,*) nyou !要素数の読み込み do i=1,nyou !各要素の頂点数?、各頂点の節点番号、カラーコード?の読み込み read(9,*) kazu,mesh(i,1),mesh(i,2),mesh(i,3),mozi end do ! do i=0,nset do j=0,nset !隣り合う節点i,jの中間の節点番号の初期値をすべて-1にしておく narabi(i,j)=-1 end do end do ! !!!!!!!!!!!! read(9,*) mozi0 !1行目の文字を読み込む read(9,*) mozi0 !1行目の文字を読み込む read(9,*) nset !節点数の読み込み do i=0,nset-1 read(9,*) xu(i),yu(i),zu(i) !各節点の座標値の読み込み end do ! do i=1,nyou !各要素の頂点数?、各頂点の節点番号、カラーコード?の読み込み read(9,*) kazu,mesh2(i,1),mesh2(i,2),mesh2(i,3),mozi end do ! do i=0,nset do j=0,nset !隣り合う節点i,jの中間の節点番号の初期値をすべて-1にしておく narabi(i,j)=-1 end do end do ! ! ! print*,'solid' do i=1,nyou !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=x(mesh(i,2)+1); y2=y(mesh(i,2)+1); z2=z(mesh(i,2)+1) x3=x(mesh(i,3)+1); y3=y(mesh(i,3)+1); z3=z(mesh(i,3)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,2)+1); y1=y(mesh(i,2)+1); z1=z(mesh(i,2)+1) x2=x(mesh(i,3)+1); y2=y(mesh(i,3)+1); z2=z(mesh(i,3)+1) x3=xu(mesh(i,3)+1); y3=yu(mesh(i,3)+1); z3=zu(mesh(i,3)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,3)+1); y1=y(mesh(i,3)+1); z1=z(mesh(i,3)+1) x2=xu(mesh(i,3)+1); y2=yu(mesh(i,3)+1); z2=zu(mesh(i,3)+1) x3=x(mesh(i,2)+1); y3=y(mesh(i,2)+1); z3=z(mesh(i,2)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=x(mesh(i,2)+1); y2=y(mesh(i,2)+1); z2=z(mesh(i,2)+1) x3=xu(mesh(i,3)+1); y3=yu(mesh(i,3)+1); z3=zu(mesh(i,3)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=x(mesh(i,2)+1); y2=y(mesh(i,2)+1); z2=z(mesh(i,2)+1) x3=xu(mesh(i,3)+1); y3=yu(mesh(i,3)+1); z3=zu(mesh(i,3)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,2)+1); y1=y(mesh(i,2)+1); z1=z(mesh(i,2)+1) x2=xu(mesh(i,3)+1); y2=yu(mesh(i,3)+1); z2=zu(mesh(i,3)+1) x3=xu(mesh(i,2)+1); y3=yu(mesh(i,2)+1); z3=zu(mesh(i,2)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=xu(mesh(i,3)+1); y1=yu(mesh(i,3)+1); z1=zu(mesh(i,3)+1) x2=xu(mesh(i,2)+1); y2=yu(mesh(i,2)+1); z2=zu(mesh(i,2)+1) x3=x(mesh(i,1)+1); y3=y(mesh(i,1)+1); z3=z(mesh(i,1)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=x(mesh(i,2)+1); y2=y(mesh(i,2)+1); z2=z(mesh(i,2)+1) x3=xu(mesh(i,2)+1); y3=yu(mesh(i,2)+1); z3=zu(mesh(i,2)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=xu(mesh(i,2)+1); y2=yu(mesh(i,2)+1); z2=zu(mesh(i,2)+1) x3=xu(mesh(i,3)+1); y3=yu(mesh(i,3)+1); z3=zu(mesh(i,3)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=xu(mesh(i,2)+1); y1=yu(mesh(i,2)+1); z1=zu(mesh(i,2)+1) x2=xu(mesh(i,3)+1); y2=yu(mesh(i,3)+1); z2=zu(mesh(i,3)+1) x3=xu(mesh(i,1)+1); y3=yu(mesh(i,1)+1); z3=zu(mesh(i,1)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=xu(mesh(i,3)+1); y1=yu(mesh(i,3)+1); z1=zu(mesh(i,3)+1) x2=xu(mesh(i,1)+1); y2=yu(mesh(i,1)+1); z2=zu(mesh(i,1)+1) x3=x(mesh(i,1)+1); y3=y(mesh(i,1)+1); z3=z(mesh(i,1)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! x1=x(mesh(i,1)+1); y1=y(mesh(i,1)+1); z1=z(mesh(i,1)+1) x2=xu(mesh(i,2)+1); y2=yu(mesh(i,2)+1); z2=zu(mesh(i,2)+1) x3=xu(mesh(i,1)+1); y3=yu(mesh(i,1)+1); z3=zu(mesh(i,1)+1) call housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) print*,'facet normal',xn,yn,zn print*,'outer loop' print*,'vertex',x1,y1,z1 print*,'vertex',x2,y2,z2 print*,'vertex',x3,y3,z3 print*,'endloop' print*,'endfacet' !!!!!!!!!!!!!!!!!!!!!! print*, 'endlsolid' end do ! ! ! end !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! subroutine daiya(r,h,m,n) implicit real*8 (a-h, o-z) pi=2.d0*asin(1.d0) nset=m*(n+1) !節点数 nyou=2*m*n !要素数 open(7,file='yousosuu.txt') write(7,*) nyou write(7,*) m,n close(7) th=2.d0*pi/m !ダイヤ1つぶんの中心角 dh=h/n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない) ! open(8,file='daiya3d.obj') write(8,'("3DG1")') write(8,*) nset do j=0,n !高さ方向に do i=0,m-1 !円周方向に1回りずつ !xyz座標を出力 write(8,*) r*cos(i*th+j*th/2.),r*sin(i*th+j*th/2.),j*dh end do end do ! !1周ずつ逆三角形と順三角形の頂点の節点を交互に外から見て左回りに出力 do j=0,n-1 !高さ方向 do i=0,m-1 !周方向に if(i/=m-1) then !周の開始から一つ前まで !逆三角形の右上頂点、左上頂点、下頂点 write(8,*) 3,j*m+i+1,j*m+i,(j+1)*m+i,' ','0xcccccc' !三角形の左下頂点、右下頂点、上 write(8,*) 3,(j+1)*m+i,(j+1)*m+i+1,j*m+i+1,' ','0xcccccc' else !周の一番最後は、開始点の頂点とくっつくように write(8,*) 3,j*m,j*m+i,(j+1)*m+i,' ','0xcccccc' write(8,*) 3,(j+1)*m+i,(j+1)*m,j*m,' ','0xcccccc' end if end do end do write(8,'("#")') !daiya.objファイルから要素数を読み取るとき用 ! ! return end ! subroutine housen(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn) implicit real*8 (a-h, o-z) ax=x2-x1 ay=y2-y1 az=z2-z1 bx=x3-x1 by=y3-y1 bz=z3-z1 ab=sqrt(ax**2+ay**2+az**2)*sqrt(bx**2+by**2+bz**2) anaib=ax*bx+ay*by+az*bz costh=anaib/ab sinth=sqrt(1.d0-costh**2) ookisa=ab*sinth xn=ax*bz-az*by yn=az*bx-ax*bz zn=ax*by-ay*bx xn=xn/ookisa yn=yn/ookisa zn=zn/ookisa return end