!saikasitaを出力ver
!片持ち梁の端を固定。
!自由端の上部に線載荷して、一番したのところのたわみの値を出す。
!14/5/16 滝田さんのccxc3d4unv.f90を少しだけ改良ver2
!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='koumoku.unv')
open(8,file='koumoku.inp')
!高さ（Salomeでのy方向の幅)(m単位で設定）
takasa=5.0d-2
!
p=1.0d-4!荷重(MN)
ell=2.0d-1!inpファイル内での梁の長さ(載荷節点用)
nsettensuu=121!節点数を入力(Node)
nedgeyouso=56!1D(edges)要素の数の入力
nfaceyouso=188!2D(faces)要素の数の入力
nvolumeyouso=404!3D(volume)要素の数の入力
!ヤング率の入力
ezz=7.694d3!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
write(8,'("** 要素番号とその要素を囲む節点番号")')
write(8,'("*ELEMENT, TYPE=C3D4, ELSET=Eall")')
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
read(7,*) nset1,nset2,nset3,nset4
write(8,*) nyou-nfaceyouso-nedgeyouso,',',nset1,',',nset2,',',nset3,',',nset4
end do
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
close(7)
!write(8,'("** 要素グループここまで")')
!write(8,'("** 拘束節点を以下に挿入")')
write(8,'("*NSET,NSET=kousoku")')
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,'("** 拘束節点ここまで")')
write(8,'("** 載荷節点を以下に挿入")')
write(8,'("*NSET,NSET=saika")')
nsaikasettensuu=0
do i=1,nsettensuu    !節点数
if (ell-0.00002<z(i)/1.0d3.and.z(i)/1.0d3<ell+0.00002&
    &.and.-0.0001<y(i).and.y(i)<0.0001) then!片持ち梁の端っこのz座標を条件にする
write(8,*) i,','
nsaikasettensuu=nsaikasettensuu+1
end if
end do
write(8,'("*NSET,NSET=saikasita")')
do i=1,nsettensuu
if (ell-0.00001<z(i)/1.0d3.and.z(i)/1.0d3<ell+0.00001d0&
&.and.y(i)/1.0d3<=takasa+0.0001.and.y(i)/1.0d3>=takasa-0.0001) then
write(8,*) i,','
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=E1,MATERIAL=EL1")')
write(8,'("*SOLID SECTION,ELSET=E2,MATERIAL=EL2")')
write(8,'("*STEP")')
write(8,'("*STATIC,SOLVER=SPOOLES")')
write(8,'("*BOUNDARY")')
write(8,'("kousoku,1,3")')
write(8,'("*CLOAD")')
! "saika"で定義した節点の2(y)方向にp/載荷節点数を載荷
write(8,'(A,a,i6,a,1pe13.5)')'saika',',',2,',',p/nsaikasettensuu
write(8,'("*NODE PRINT,NSET=saikasita")')
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)
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
