! EUC-jp
! 文字コードはEUC-jp
! 自由形式のfortranで書いています。
! g77でコンパイルする場合は、
!g77 -ffree-form -o smv smv.f
! のようにコンパイルしてください。
! gfortran でコンパイルする場合は文字コードをutf-8にしてください。
!一枚の板の大きさに統一した場合の
!ダイヤカット円筒の線形座屈解析を有限要素法で解くプログラム
! 09/3/22版
! プログラムを主に作った人：後藤文彦，工藤康広
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! このプログラムは無保証です。
! このプログラムは改造／再配布して構いません。
! 但し、私的利用に留まらない場合は、もとのプログラムが
! http://www.str.ce.akita-u.ac.jp/~gotou/programoj/
! から取ってきたものと分かるように明示して下さい。
! 改造／再配布後も上記の条件を保持して下さい。
! バグの報告は、
! http://www.str.ce.akita-u.ac.jp/~gotou/png/meeru.png
! にお願いします。
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!program daiya
implicit real*8 (a-h, o-z)
pi=2.d0*asin(1.d0)
!!半径、高さ、周方向分割数、高さ方向の三角形の個数を設定して
hankei=3.d-2 !半径
yoko=2*hankei*pi !円周
tate=1.2d-1 !高さ 
m=8!求めたい周方向パターン数を2倍して書き込む(8って書いたらパターン数は4になってる)
!read*,n3kaku
n3kaku=4!高さ方向パターン数
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!ここで荷重を入力(荷重を載荷節点数で割る)
saika=1.d-4/m/(n3kaku+1)
!載荷方向を設定します
nhoukou=-1 !(中心方向に圧縮載荷したい場合)
!nhoukou=1 (中心に引張載荷したい場合)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
sykakudo=2*pi/m
r=yoko/m/sin(sykakudo)
syahen=(tate/n3kaku)
!syahen=(tate/n3kaku)
tobi=r*(1-cos(sykakudo))
!tobi=r*(1-cos(sykakudo))
h=sqrt(syahen**2-tobi**2)
!h=sqrt(-1*(syahen**2-tobi**2))
!tamesi=sqrt(h**2+tobi**2)
!print*, 'r',r
!print*, 'sya',syahen
!print*, 'tobi',tobi
!print*, 'h',h
!print*, 'tame',tamesi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
n=n3kaku*2 !高さ方向の分割数(一つの三角形の高さ方向は二分割してるから)
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/2 !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない)
!
!
!
open(8,file='syutkotei.obj')
write(8,'("3DG1")')
write(8,*) nset
!
!
do j=0,n-2,2
!
!
!一段目
ensyuu=1
do i=0,m-1 !円周方向に1回りずつ
!xyz座標を出力
if(ensyuu>0) then
write(8,*) r*cos(th)*cos(i*th+j*th/2),&
&r*cos(th)*sin(i*th+j*th/2),j*dh
else
write(8,*) r*cos(i*th+j*th/2),r*sin(i*th+j*th/2),j*dh
end if
ensyuu=ensyuu*(-1.)
end do
!
!
!
!二段目
ensyuu=1
do i=0,m-1 !円周方向に1回りずつ
if(ensyuu>0) then
write(8,*) (r*cos((i+1)*th+j*th/2)+r*cos((i-1)*th+(j+2)*th/2.))/2,&
&(r*sin((i+1)*th+j*th/2.)+r*sin((i-1)*th+(j+2)*th/2.))/2,(j+1)*dh
else
write(8,*) (r*cos(i*th+j*th/2)+r*cos(i*th+(j+2)*th/2))/2,&
&(r*sin(i*th+j*th/2)+r*sin(i*th+(j+2)*th/2))/2,(j+1)*dh
end if
ensyuu=ensyuu*(-1.)
end do
end do
!
!最終段
ensyuu=1
do i=0,m-1 !円周方向に1回りずつ
!xyz座標を出力
if(ensyuu>0) then
write(8,*) r*cos(th)*cos(i*th+n*th/2),&
&r*cos(th)*sin(i*th+n*th/2),n*dh
else
write(8,*) r*cos(i*th+n*th/2),r*sin(i*th+n*th/2),n*dh
end if
ensyuu=ensyuu*(-1.)
end do
!
!
!
!1周ずつ逆三角形と順三角形の頂点の節点を交互に外から見て左回りに出力
!do j=0,n-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ファイルから要素数を読み取るとき用
close(8)
!
!
open(9,file='yokoka.inp')
write(9,*) m*(n3kaku+1)*2
do j=0,n-2,2
!一段目
ensyuu=1
do i=0,m-1 !円周方向に1回りずつ
!xy座標を出力
write(9,*) m*j+i+1,',',1,',',nhoukou*saika*cos(i*th+j*th/2)
write(9,*) m*j+i+1,',',2,',',nhoukou*saika*sin(i*th+j*th/2)
end do
end do
!
!最終段
!ensyuu=1
do i=0,m-1 !円周方向に1回りずつ
write(9,*) m*j+i+1,',',1,',',nhoukou*saika*cos(i*th+n*th/2)
write(9,*) m*j+i+1,',',2,',',nhoukou*saika*sin(i*th+n*th/2)
end do
close(9)
!
end
