!竜太さん用
!Salomeで円柱の円の中心を原点にして作ったものに特化
!c3d4,c3d6,c3d8が混在しててもok
!グルーピングには非対応,したい人は各自で書いて
!c3d6,c3d4混在用を改良して、c3d8が混ざってても良くなった
!salomeではpyramidができているとccxで解けなくなるので注意!
!Salome6.6用にSalome内でmm単位で書くと、inpファイルでm単位
!下でヤング率とポアソン比設定
implicit real*8(a-h,o-z)
dimension x(10000000),y(1000000),z(1000000),e66(6,6),b(3,3)
open(7,file='en15.unv')
open(8,file='en15.inp')
!一時ファイルを初期化
open(9,file='xpmvdc3d6.inp', status='replace')
open(10,file='xpmvdc3d8.inp', status='replace')
close(9)
close(10)
!以下はSalomeでメッシュを切った(compute)した時にでる
!値の入力(ツリーのMesh_1を右クリック,MeshInformationでも見れる)
!Mesh Information(0の時は0って入れて)
nsettensuu=20449!節点数を入力(Node)
nedgeyouso=358!1D(edges)要素の数の入力
nfaceyouso=26386!2D(faces)要素の数の入力
nvolumeyouso=83661!3D(volume)要素の数の入力
nhex=0!hexahedrons(C3D8要素)
nprism=0!prism(C3D6要素)
!載荷拘束用##############################
!Salomeでのy方向の幅(m単位で設定）
enr=1.6d-1 !(半径)
enr2=enr*2 !(直径)
p=1.0d-3!荷重(MN)
!inpファイル内での梁の長さ(Salomeでのz方向の長さ)(m単位)
ell=2.0d0
!以下直行異方性のための入力################33
!ヤング率の入力
ezz=5.00d3!z軸方向のヤング率
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!#############################################
write(8,'("** 節点番号と座標")')
write(8,'("*NODE, NSET=Nall")')
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
read (7,'()') 
do i=1,nsettensuu    !節点数
read (7,*) nset,n1,n2,n3
read(7,*) x(i),y(i),z(i)
write(8,'(i6,a,1pe13.5,a,1pe13.5a,1pe13.5)') nset,',',x(i)/1.0d3&
&,',',y(i)/1.0d3,',',z(i)/1.0d3
end do
!ここまでok
write(8,'("** 要素番号とその要素を囲む節点番号")')
write(8,'("*ELEMENT, TYPE=C3D4, ELSET=E1")')
read (7,'()') 
read (7,'()') 
read (7,'()') 
do i=1,nedgeyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
do i=1,nfaceyouso
read(7,*) nyou,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
end do
do i=1,nvolumeyouso
read(7,*) nyou,n1,n2,n3,n4,n5
if (n5 == 4) then!C3D4の時
read(7,*) nset1,nset2,nset3,nset4
write(8,*) nyou-nfaceyouso-nedgeyouso,',',nset1,',',nset2,',',nset3,',',nset4
else if (n5 == 6) then!C3D6の時
open(9,file='xpmvdc3d6.inp')
read(7,*) nset1,nset2,nset3,nset4,nset5,nset6
write(9,*) nyou-nfaceyouso-nedgeyouso,',',nset1,',',nset2,',',nset3,',',nset4&
    &,',',nset5,',',nset6
else if (n5 == 8) then!C3D8の時
open(10,file='xpmvdc3d8.inp')
read(7,*) nset1,nset2,nset3,nset4,nset5,nset6,nset7,nset8
write(10,*) nyou-nfaceyouso-nedgeyouso,',',nset1,',',nset2,',',nset3,',',nset4&
    &,',',nset5,',',nset6,',',nset7,',',nset8
end if
end do
close(7)
close(9)
close(10)
if (nprism>1) then
open(9,file='xpmvdc3d6.inp')
write(8,'("*ELEMENT, TYPE=C3D6, ELSET=Ec3d6")')
do i=1,nprism
read(9,*) nyoup,n1,n2,n3,n4,n5,n6
write(8,*) nyoup,',',n1,',',n2,',',n3,',',n4&
    &,',',n5,',',n6
end do
end if
if (nhex>1) then
open(10,file='xpmvdc3d8.inp')
write(8,'("*ELEMENT, TYPE=C3D8, ELSET=Ec3d8")')
do i=1,nhex
read(10,*) nyoup,n1,n2,n3,n4,n5,n6,n7,n8
write(8,*) nyoup,',',n1,',',n2,',',n3,',',n4&
    &,',',n5,',',n6,',',n7,',',n8
end do
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!write(8,'("**ここから1個目の要素グループをよんでかく")')
!write(8,'("** 要素グループを以下に挿入")')
!write(8,'("*ELSET,ELSET=E1")')
!read (7,'()')
!read (7,'()')
!read (7,'()')
!read(7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!!要素数が偶数
!if (mod (nyou8,2) == 0) then
!read (7,'()')
!do i=1,nyou8/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!write(8,*) n6-nfaceyouso-nedgeyouso,','
!end do
!!要素数が奇数の場合
!else
!read (7,'()')
!do i=1,(nyou8-1)/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!write(8,*) n6-nfaceyouso-nedgeyouso,','
!end do
!read(7,*) n1,n2,n3,n4
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!end if
!write(8,'("**ここから２個目の要素グループをよんでかく")')
!write(8,'("*ELSET,ELSET=E2")')
!read (7,*) nyou1,nyou2,nyou3,nyou4,nyou5,nyou6,nyou7,nyou8
!if (mod (nyou8,2) == 0)   then!要素数が偶数の場合
!read (7,'()')
!do i=1,nyou8/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!write(8,*) n6-nfaceyouso-nedgeyouso,','
!end do
!else!要素数が奇数の場合
!read (7,'()')
!do i=1,(nyou8-1)/2
!read(7,*) n1,n2,n3,n4,n5,n6,n7,n8
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!write(8,*) n6-nfaceyouso-nedgeyouso,','
!end do
!read(7,*) n1,n2,n3,n4
!write(8,*) n2-nfaceyouso-nedgeyouso,','
!end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(8,'("** 要素グループここまで")')
write(8,'("** 拘束節点を以下に挿入")')
write(8,'("*NSET,NSET=kousoku1")')!一番した全部xyz全部拘束
do i=1,nsettensuu    !節点数
if (-0.0001<z(i).and.z(i)<0.0001) then!ここはz(i)にd3をつけなくても良い
write(8,*) i,','
end if
end do
write(8,'("*NSET,NSET=kousoku2")')!一番上全部をxとzだけ拘束
do i=1,nsettensuu    !節点数
if (ell-0.00002<z(i)/1.0d3.and.z(i)/1.0d3<ell+0.00002)then
write(8,*) i,','
end if
end do
write(8,'("*NSET,NSET=kousoku2")')!一番上全部をxとzだけ拘束
do i=1,nsettensuu    !節点数
if (ell-0.00002<z(i)/1.0d3.and.z(i)/1.0d3<ell+0.00002)then
write(8,*) i,','
end if
end do
write(8,'("** 拘束節点ここまで")')
write(8,'("** 載荷節点を以下に挿入")')
write(8,'("*NSET,NSET=saika")')
nsaikasettensuu=0
do i=1,nsettensuu    !節点数
if (0<y(i)/1.0d3.and.y(i)/1.0d3<enr+0.0000001)then
write(8,*) i,','
nsaikasettensuu=nsaikasettensuu+1
end if
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(8,'("*MATERIAL,NAME=EL1")')
write(8,'("*ELASTIC,TYPE=ORTHOROPIC")')
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
    & gxz,gxy,poiyx,poizx,poizy)
write(8,'(f12.4,a,f12.4,a,f12.4,a,&
    &f12.4,a,f12.4,a,f12.4,a,f12.4,&
    & a ,f12.4,a)')&
    &e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
    &,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
write(8,'(f12.4,a,a)')e66(4,4),',','293'
!write(8,'("*MATERIAL,NAME=EL2")')
!write(8,'("*ELASTIC")')
!write(8,'("206.0d3,0.3")')
!write(8,'("*SOLID SECTION,ELSET=Ec3d8,MATERIAL=EL1")')
write(8,'("*SOLID SECTION,ELSET=E1,MATERIAL=EL1")')
!write(8,'("*SOLID SECTION,ELSET=Ec3d6,MATERIAL=EL2")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku1,1,3")')
write(8,'("kousoku2,1")')
write(8,'("kousoku2,3")')
write(8,'("*CLOAD")')
! "saika"で定義した節点の3(z)方向にp/載荷節点数を載荷
write(8,'(A,a,i6,a,1pe13.5)')'saika',',',2,',',p/nsaikasettensuu
write(8,'("*NODE PRINT,NSET=Nall")')
write(8,'("U")')
write(8,'("*EL PRINT,ELSET=Eall")')
write(8,'("S")')
write(8,'("*NODE FILE")')
write(8,'("U")')
write(8,'("*EL FILE,POSITION=AVERAGED AT NODES")')
write(8,'("S")')
write(8,'("*END STEP")')
close(8)
close(9)
end

!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
        &gxz,gxy,poiyx,poizx,poizy)
    implicit real*8(a-h,o-z)
    dimension e66(6,6),b(3,3)
    !応力－ひずみ行列：e66
    do i=1,6
    do j=1,6
    e66(i,j)=0.d0
    end do
    end do
    !  
    e66(1,1)=1.d0/exx
    e66(1,2)=-poixy/exx
    e66(1,3)=-poixz/exx
    e66(2,1)=-poiyx/eyy
    e66(2,2)=1.d0/eyy
    e66(2,3)=-poiyz/eyy
    e66(3,1)=-poizx/ezz
    e66(3,2)=-poizy/ezz
    e66(3,3)=1.d0/ezz
    !
    !****** 逆行列を取って、応力－ひずみ行列に変換する******
    det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
        &+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
        &-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
    b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
    b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
    b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
    b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
    b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
    b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
    b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
    b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
    b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
    do i=1,3
    do j=1,3
    e66(i,j)=b(j,i)/det
    end do
    end do
    !
    !****** せん断補正係数を入力 ******
    !print*,'xy面のせん断弾性係数を入れて下さい'
    !read*, gxy
    e66(4,4)=gyz
    !print*,'xz面のせん断弾性係数を入れて下さい'
    !read*, gxz
    e66(5,5)=gxz
    !print*,'yz面のせん断弾性係数を入れて下さい'
    !read*, gyz
    e66(6,6)=gxy
    !
    !
    return
    end
