!片持ち梁のねじり(CCX)
! 円柱(C3D6+C3D8)、弾性体、等方性材料

implicit real*8(a-h,o-z)
dimension x(1000000),y(1000000),z(1000000),e66(6,6),b(3,3)&
&,mesh(30000000,8)
!プログラムはメートル、MPa
!寸法
d=13.d-3!直径(m)
ell=3.0d0!長さ(m)
to=1.0d-6!荷重 モーメント力(MN・m)
p=0.0d-6!荷重 先端の引張力(MN)
!rcz方向の分割数
nr=4!半径方向
nc=16!円周方向
nz=2000!軸方向
!1つあたりの要素の長さとか角度とか
pi=asin(1.d0)*2.d0
th=pi*2.d0/nc!角度
rn=d/(2*nr)!半径方向
zn=ell/nz!軸方向
!材料定数
e=206.d3!ヤング率(MN)
nu=0.3!ポアソン比
g=e/(2.d0*(1+nu))!せん断弾性係数(MN)
cj=pi/2*(d/2)**4.!ねじり定数
!ねじり角
phi=to/(g*cj)*ell
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
nset=0
do k=0,nz
nset=nset+1
x(nset)=0.d0
y(nset)=0.d0
z(nset)=k*zn
do j=1,nr
do i=0,nc-1
nset=nset+1
x(nset)=rn*j*cos(th*i)
y(nset)=rn*j*sin(th*i)
z(nset)=k*zn
end do
end do
end do 
nsettenzenbu=nset!全節点数
nitimen=1+nr*nc!断面の節点数
!i5~6が整数,aが文字,f9,6が小数.1peとかはEの＊乗
!接点番号と座標の出力
print'("*NODE, NSET=Nall")'
do i=1,nsettenzenbu
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
    &i,',',x(i),',',y(i),',',z(i)
end do
i=nsettenzenbu+1
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
    &i,',',0.d0,',',0.d0,',',0.d0
i=nsettenzenbu+2
print'(i6,a,1pe17.9,a,1pe17.9,a,1pe17.9)',&
!    &i,',',d/2.d0,',',d/2.d0,',',ell
    &i,',',0.d0,',',0.d0,',',ell
!要素番号とその要素の6接点の節点番号
nyou=0
do k=0,nz-1
do i=0,nc-2
nyou=nyou+1
mesh(nyou,1)=1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+i+1
mesh(nyou,3)=mesh(nyou,2)+1
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,2)+nitimen
mesh(nyou,6)=mesh(nyou,3)+nitimen
end do
i=nc-1
nyou=nyou+1
mesh(nyou,1)=1+k*nitimen
mesh(nyou,2)=mesh(nyou,1)+i+1
mesh(nyou,3)=mesh(nyou,2)+1-nc
mesh(nyou,4)=mesh(nyou,1)+nitimen
mesh(nyou,5)=mesh(nyou,2)+nitimen
mesh(nyou,6)=mesh(nyou,3)+nitimen
end do
nyounaka=nyou
!要素の出力
print'("*ELEMENT, TYPE=C3D6, ELSET=Enaka")'
do i=1,nyounaka
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &i,',',&
    &mesh(i,1),',',&
    &mesh(i,2),',',&
    &mesh(i,3),',',&
    &mesh(i,4),',',&
    &mesh(i,5),',',&
    &mesh(i,6)
end do
!要素番号とその要素の8接点の節点番号
do k=0,nz-1
do j=0,nr-2
do i=1,nc-1
nyou=nyou+1
mesh(nyou,1)=1+nitimen*k+nc*j+i
mesh(nyou,2)=mesh(nyou,1)+nc
mesh(nyou,3)=mesh(nyou,2)+1
mesh(nyou,4)=mesh(nyou,3)-nc
mesh(nyou,5)=mesh(nyou,1)+nitimen
mesh(nyou,6)=mesh(nyou,2)+nitimen
mesh(nyou,7)=mesh(nyou,3)+nitimen
mesh(nyou,8)=mesh(nyou,4)+nitimen
end do
i=nc
nyou=nyou+1
mesh(nyou,1)=1+nitimen*k+nc*j+i
mesh(nyou,2)=mesh(nyou,1)+nc
mesh(nyou,3)=mesh(nyou,2)+1-nc
mesh(nyou,4)=mesh(nyou,3)-nc
mesh(nyou,5)=mesh(nyou,1)+nitimen
mesh(nyou,6)=mesh(nyou,2)+nitimen
mesh(nyou,7)=mesh(nyou,3)+nitimen
mesh(nyou,8)=mesh(nyou,4)+nitimen
end do
end do
!要素の出力
print'("*ELEMENT, TYPE=C3D8, ELSET=Esoto")'
do i=nyounaka+1,nyou
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
    &i,',',&
    &mesh(i,1),',',&
    &mesh(i,2),',',&
    &mesh(i,3),',',&
    &mesh(i,4),',',&
    &mesh(i,5),',',&
    &mesh(i,6),',',&
    &mesh(i,7),',',&
    &mesh(i,8)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!固定端の節点群を書きだす
print'("*NSET,NSET=kotei")'
do i=1,nitimen
nkotei=i
print'(i6,a)',nkotei
end do
!載荷点の節点群を書きだす
print'("*NSET,NSET=saika")'
do i=1,nitimen
nsaika=nitimen*nz+i
print'(i6,a)',nsaika
end do
!載荷点の回転軸を決める
print'(A,i6,a,A,i6)',&
    &'*RIGID BODY,NSET=saika,REF NODE=',&
    &nsettenzenbu+2,',',&
    &'ROT NODE=',&
    &nsettenzenbu+1
!境界条件
print'("*BOUNDARY")'
print'("kotei,1,3")'
print'(i6,a,i6,a,i6,a,1pe13.5)',nsettenzenbu+2,',',1,',',1,',',0.0
print'(i6,a,i6,a,i6,a,1pe13.5)',nsettenzenbu+2,',',2,',',2,',',0.0
! 木材の設定
print'("*MATERIAL,NAME=EL1")'
print'("*ELASTIC")'
print'("206d3, 0.3d0")'
print'("*SOLID SECTION,ELSET=Enaka,MATERIAL=EL1")'
print'("*SOLID SECTION,ELSET=Esoto,MATERIAL=EL1")'
! 載荷条件
!print'("*STEP,NLGEOM,INC=1000")'
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
!print'("0.1,1.")'
print'("*BOUNDARY")'
print'(i6,a,i6,a,i6,a,1pe13.5)',nsettenzenbu+1,',',3,',',3,',',phi
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',3,',',p/nitimen
print'("*BUCKLE")'
print'("10,0.01")'
! datのオプション
print'("*NODE PRINT,NSET=Nall")'
print'("U")'
!print'("*EL PRINT,ELSET=Enaka")'
!print'("S")'
!print'("*EL PRINT,ELSET=Esoto")'
!print'("S")'
print'("*NODE FILE")'
print'("U")'
print'("*EL FILE")'
print'("S")'
print'("*END STEP")'
!
end

