! 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 

