!moku3d10.f90を修正して、接触解析にも対応させてみる。
!斜めになった部材とかにも木材の異方性を設定できる。

! 拘束とか載荷の節点群は、座標指定ではなく、Salomeのmeshのグループで指定する。
! なので、Salomeで、Mesh_1を右クリック-->Create Groups from Geometry-->Nodes
! から、geometryで予めcreate groupしたグループを選択して作っておく。
! 材料別の要素群も同様に、meshのグループを作る必要があるが、
! 要素は並び順情報が要るので、Create Groups from Geometry-->Elements で選択する。

! ちなみに、  /#  で 修正必要箇所 に跳べる。

implicit real*8(a-h,o-z)
dimension x(500000),y(500000),z(500000)
dimension e66(6,6),b(3,3),mesh(5000000,11),men(100000,6)
dimension t(3,3),alp(4),bt(4),gmm(4)

!---各種入力データ---
!#--入出力ファイル名
open(7,file='ataru.unv')
open(8,file='ataru.inp')
!#--荷重
p=1.d-6!荷重(MN)
!#--Mesh_1 の Mesh Information を見て入力
nnode=296   !2次要素の節点数(Nodes)
nedge=32    !1D(edges)の数
nface=76    !2D(faces)の数
nvolume=148 !3D(volumes)の数
!---sub mesh の Mesh Information を見て入力
!#--各部材の要素数(volumes)と固定軸回りのオイラー角？
nbuzai=1 !木材の向きが何通りあるか
pi=asin(1.d0)*2.d0
nyouso1=136 !木材の要素数 
alp(1)= 1.d0/2*pi     ; bt(1)= 0.d0         ; gmm(1)=0.d0
!nyouso2=14790
!alp(2)=11.d0/6.d0*pi ; bt(2)= 0.d0         ; gmm(2)=0.d0        
!nyouso3=14655
!alp(3)= 1.d0/6.d0*pi ; bt(3)= 0.d0         ; gmm(3)=0.d0        
!nyouso4=22244
!alp(4)= 3.d0/2.d0*pi ; bt(4)= 0.d0         ; gmm(4)=0.d0        
nyouso0=12    !鋼材
!#--境界条件の節点
nkotei1=41     !固定グループの節点数(Nodes)
!nkotei2=12     !固定グループの節点数(Nodes)
nsaika=9     !載荷グループの節点数(Nodes)
!#--接触面の面数
nslave=2     !スレーブ面の面数
nmaster=16   !マスター面の面数
!#--変位を読み取る節点(あると便利）
!nyomi1=8       !読み取りグループの節点数(Nodes)
!nyomi2=7       !読み取りグループの節点数(Nodes)
!#--材料定数
!---鋼材と木材の摩擦係数
cf=0.4
!---鋼材のヤング率(MPa)
young=206.d3
!young=2.84d3
!---鋼材のポアソン比
poisson=0.3d0
!---木材のヤング率(MPa)
ezz=7.5d3
exx=ezz/25.0d0
eyy=ezz/25.0d0
!---木材のポアソン比
poizx=0.4d0
poixy=poizx/ezz*exx
poixz=poixy
poiyx=poizx/ezz*eyy
poiyz=poiyx
!---せん断弾性係数(MPa)
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!---

!---inpファイルの節点情報(*NODE)を作る---
write(8,'("*NODE, NSET=Nall")')
do i=1,19
read (7,'()')!ヘッダーの読み飛ばし
end do
!---節点番号と節点座標の読み込みと書き出し
do i=1,nnode
read (7,*) nset,n1,n2,n3
read(7,*) x(i),y(i),z(i)
write(8,'(i10,a,1pe17.9,a,1pe17.9a,1pe17.9)') nset,',',x(i)/1.0d3&
&,',',y(i)/1.0d3,',',z(i)/1.0d3
end do
!---

!---inpファイルの要素情報を読み取る---
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---Edges情報の読み飛ばし
do i=1,nedge
read(7,*) nset,n1,n2,n3,n4,n5
read(7,*) n6,n7,n8
read(7,*) n9,n10
end do
!---Faces情報の読み取り
do i=1,nface
read(7,*) nset,n1,n2,n3,n4,n5
!read(7,*) men(nset,1),men(nset,3),n5,men(nset,2),n4,n6
read(7,*) men(nset,1),n2,men(nset,2),n5,men(nset,3),n6
men(nset,4)=men(nset,1)
men(nset,5)=men(nset,2)
men(nset,6)=men(nset,3)
end do
!---Volumes情報の読み込みと書き出し、ついでに体積計算
taiseki=0.d0
do i=1,nvolume
read(7,*) nset,n1,n2,n3,n4,n5
mesh(i,11)=nset !unvの要素番号を記憶
read(7,*) mesh(nset,1),mesh(nset,5),mesh(nset,2),mesh(nset,6)&
    &,mesh(nset,3),mesh(nset,7),mesh(nset,8),mesh(nset,9)
read(7,*) mesh(nset,10),mesh(nset,4)
!---体積の計算
x1=x(mesh(nset,1));y1=y(mesh(nset,1));z1=z(mesh(nset,1))
x2=x(mesh(nset,2));y2=y(mesh(nset,2));z2=z(mesh(nset,2))
x3=x(mesh(nset,3));y3=y(mesh(nset,3));z3=z(mesh(nset,3))
x4=x(mesh(nset,4));y4=y(mesh(nset,4));z4=z(mesh(nset,4))
taiseki=taiseki+(&
& (x4-x1)*((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)) &
&+(y4-y1)*((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)) &
&+(z4-z1)*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)) &
&)/6
end do
print*,'体積は',taiseki
read (7,'()') 
read (7,'()') 
read (7,'()') 
!---

!---inpファイルの節点群(*NSET)を作る---
!   以下のグループの並び順は、Salomeで sub mesh を作った順(unvの順)に並び替える
!#--載荷するグループの節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=saika")')
do i=1,int(nsaika/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nsaika,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!#--拘束するグループ1の節点番号の読み込みと書き出し
read (7,'()') 
read (7,'()') 
write(8,'("*NSET,NSET=kousoku1")')
do i=1,int(nkotei1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
write(8,*) nset1,','
write(8,*) nset2,','
enddo
if (mod(nkotei1,2)==1) then
read(7,*) n1,nset1,n2,n3
write(8,*) nset1,','
endif
!!#--拘束するグループ2の節点番号の読み込みと書き出し
!read (7,'()') 
!read (7,'()') 
!write(8,'("*NSET,NSET=kousoku2")')
!do i=1,int(nkotei2/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!write(8,*) nset1,','
!write(8,*) nset2,','
!enddo
!if (mod(nkotei2,2)==1) then
!read(7,*) n1,nset1,n2,n3
!write(8,*) nset1,','
!endif
!!   読み取りグループは、あると便利
!!#--読み取りグループ1の節点番号の読み込みと書き出し
!read (7,'()') 
!read (7,'()') 
!write(8,'("*NSET,NSET=yomi1")')
!do i=1,int(nyomi1/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!write(8,*) nset1,','
!write(8,*) nset2,','
!enddo
!if (mod(nyomi1,2)==1) then
!read(7,*) n1,nset1,n2,n3
!write(8,*) nset1,','
!endif
!!#--読み取りグループ2の節点番号の読み込みと書き出し
!read (7,'()') 
!read (7,'()') 
!write(8,'("*NSET,NSET=yomi2")')
!do i=1,int(nyomi2/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!write(8,*) nset1,','
!write(8,*) nset2,','
!enddo
!if (mod(nyomi2,2)==1) then
!read(7,*) n1,nset1,n2,n3
!write(8,*) nset1,','
!endif
!---

!---inpファイルの要素群(*ELSET)を作る---
nyou=0
!#--材料1
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E1")')
do i=1,int(nyouso1/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso1,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!!#--材料2
!read (7,'()') 
!read (7,'()') 
!write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E2")')
!do i=1,int(nyouso2/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
!    &,mesh(nset2,3),',',mesh(nset2,4),','&
!    &,mesh(nset2,5),',',mesh(nset2,6),','&
!    &,mesh(nset2,7),',',mesh(nset2,8),','&
!    &,mesh(nset2,9),',',mesh(nset2,10),','
!enddo
!if (mod(nyouso2,2)==1) then
!read(7,*) n1,nset1,n2,n3
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!endif
!!#--材料3
!read (7,'()') 
!read (7,'()') 
!write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E3")')
!do i=1,int(nyouso3/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
!    &,mesh(nset2,3),',',mesh(nset2,4),','&
!    &,mesh(nset2,5),',',mesh(nset2,6),','&
!    &,mesh(nset2,7),',',mesh(nset2,8),','&
!    &,mesh(nset2,9),',',mesh(nset2,10),','
!enddo
!if (mod(nyouso3,2)==1) then
!read(7,*) n1,nset1,n2,n3
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!endif
!!#--材料4
!read (7,'()') 
!read (7,'()') 
!write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E4")')
!do i=1,int(nyouso4/2)
!read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
!    &,mesh(nset2,3),',',mesh(nset2,4),','&
!    &,mesh(nset2,5),',',mesh(nset2,6),','&
!    &,mesh(nset2,7),',',mesh(nset2,8),','&
!    &,mesh(nset2,9),',',mesh(nset2,10),','
!enddo
!if (mod(nyouso4,2)==1) then
!read(7,*) n1,nset1,n2,n3
!nyou=nyou+1
!write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
!    &,mesh(nset1,3),',',mesh(nset1,4),','&
!    &,mesh(nset1,5),',',mesh(nset1,6),','&
!    &,mesh(nset1,7),',',mesh(nset1,8),','&
!    &,mesh(nset1,9),',',mesh(nset1,10),','
!endif
!#--材料0
read (7,'()') 
read (7,'()') 
write(8,'("*ELEMENT, TYPE=C3D10, ELSET=E0")')
do i=1,int(nyouso0/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
nyou=nyou+1
write(8,*) nyou,',',mesh(nset2,1),',',mesh(nset2,2),','&
    &,mesh(nset2,3),',',mesh(nset2,4),','&
    &,mesh(nset2,5),',',mesh(nset2,6),','&
    &,mesh(nset2,7),',',mesh(nset2,8),','&
    &,mesh(nset2,9),',',mesh(nset2,10),','
enddo
if (mod(nyouso0,2)==1) then
read(7,*) n1,nset1,n2,n3
nyou=nyou+1
write(8,*) nyou,',',mesh(nset1,1),',',mesh(nset1,2),','&
    &,mesh(nset1,3),',',mesh(nset1,4),','&
    &,mesh(nset1,5),',',mesh(nset1,6),','&
    &,mesh(nset1,7),',',mesh(nset1,8),','&
    &,mesh(nset1,9),',',mesh(nset1,10),','
endif
!---接触面
!#--マスター面
read (7,'()') 
read (7,'()') 
write(8,'("*SURFACE,NAME=SM")')
do i=1,int(nmaster/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
do j=1,nvolume
do k=1,3
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset1,k)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset2,k)==mesh(mesh(j,11),2).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
enddo
enddo
enddo
if (mod(nmaster,2)==1) then
read(7,*) n1,nset1,n2,n3
do j=1,nvolume
do k=1,3
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset1,k)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
enddo
enddo
endif
!#--スレーブ面
read (7,'()') 
read (7,'()') 
write(8,'("*SURFACE,NAME=SS")')
do i=1,int(nslave/2)
read(7,*) n1,nset1,n2,n3,n4,nset2,n5,n6
do j=1,nvolume
do k=1,3
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset1,k)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset2,k)==mesh(mesh(j,11),2).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset2,k)==mesh(mesh(j,11),1).and.&
   &men(nset2,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset2,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
enddo
enddo
enddo
if (mod(nslave,2)==1) then
read(7,*) n1,nset1,n2,n3
do j=1,nvolume
do k=1,3
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),2)) then
write(8,*) j,',',"S1"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S2"
endif
if (men(nset1,k)==mesh(mesh(j,11),2).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),3).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),4)) then
write(8,*) j,',',"S3"
endif
if (men(nset1,k)==mesh(mesh(j,11),1).and.&
   &men(nset1,k+1)==mesh(mesh(j,11),4).and.&
   &men(nset1,k+2)==mesh(mesh(j,11),3)) then
write(8,*) j,',',"S4"
endif
enddo
enddo
endif
!---
close(7)

!#--境界条件(拘束方向は、x:1,y:2,z:3)
write(8,'("*BOUNDARY")')
write(8,'("kousoku1,1,3")')
write(8,'("saika,1,2")')
!write(8,'("kousoku2,1")')
!write(8,'("kousoku2,3")')

!#--接触設定
!write(8,'("*CONTACT PAIR,INTERACTION=SI1")')!CCX2.6
write(8,'("*CONTACT PAIR,INTERACTION=SI1,TYPE=NODE TO SURFACE")')!CCX2.9
write(8,'("SS,SM")')!スレーブ面,マスター面
write(8,'("*SURFACE INTERACTION,NAME=SI1")')
write(8,'("*SURFACE BEHAVIOR,PRESSURE-OVERCLOSURE=EXPONENTIAL")')
write(8,'("1.e-10,1.e-10")')!接触と判定する距離(m),距離0の時の応力?
write(8,'("*FRICTION")')
write(8,'(f12.4,a,f12.4)')cf,',',ezz!摩擦係数,?

!#--材料の設定0
write(8,'("*MATERIAL,NAME=ELkou")')
write(8,'("*ELASTIC")')
write(8,'(1pe17.9,a,1pe17.9)')young,',',poisson
write(8,'("*SOLID SECTION,ELSET=E0,MATERIAL=ELkou")')
!   ここから、異方性の設定
!#--材料の設定:木材1
nb=1!部材番号
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
write(8,'("*MATERIAL,NAME=ELmoku")')
write(8,'("*ELASTIC,TYPE=ORTHO")')
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)'),e66(4,4)
write(8,'("*ORIENTATION,NAME=OR1,SYSTEM=RECTANGULAR")')
write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
write(8,'("*SOLID SECTION,ELSET=E1,MATERIAL=ELmoku,ORIENTATION=OR1")')
!!#--材料の設定:木材2
!nb=2!部材番号
!call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
!& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
!write(8,'("*ORIENTATION,NAME=OR2,SYSTEM=RECTANGULAR")')
!write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
!&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
!write(8,'("*SOLID SECTION,ELSET=E2,MATERIAL=Emoku,ORIENTATION=OR2")')
!!#--材料の設定:木材3
!nb=3!部材番号
!call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
!& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
!write(8,'("*ORIENTATION,NAME=OR3,SYSTEM=RECTANGULAR")')
!write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
!&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
!write(8,'("*SOLID SECTION,ELSET=E3,MATERIAL=Emoku,ORIENTATION=OR3")')
!!#--材料の設定:木材4
!nb=4!部材番号
!call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
!& gxz,gxy,poiyx,poizx,poizy,t,alp,bt,gmm,nb)
!write(8,'("*ORIENTATION,NAME=OR4,SYSTEM=RECTANGULAR")')
!write(8,'(f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9,a,f13.9)'),&
!&t(1,1),',',t(2,1),',',t(3,1),',',t(1,2),',',t(2,2),',',t(3,2)
!write(8,'("*SOLID SECTION,ELSET=E4,MATERIAL=Emoku,ORIENTATION=OR4")')
!!#--出力ステップ数
nstep=10
write(8,'("*TIME POINTS,NAME=T1")')
do i=1,nstep
ts=real(i)/nstep
write(8,'(f12.4,a)')ts,','
enddo
!---計算の設定(*STEP~*END STEP)
write(8,'("*STEP,NLGEOM")')!幾何学非線形
write(8,'("*STATIC")')
!#--荷重条件(載荷方向は、x:1,y:2,z:3)
write(8,'("*CLOAD")')
write(8,'(A,a,i6,a,1pe25.13)')'saika',',',3,',',-p/nsaika
!#--出力の設定
!write(8,'("*NODE PRINT,NSET=yomi1,TIMEPOINTS=T1")')
!write(8,'("U")')
!write(8,'("*NODE PRINT,NSET=yomi2,TIMEPOINTS=T1")')
!write(8,'("U")')
!write(8,'("*EL PRINT,GLOBAL=YES,ELSET=Eall,TIMEPOINTS=T1")')
!write(8,'("S")')
write(8,'("*NODE FILE,TIMEPOINTS=T1")')
write(8,'("U")')
write(8,'("*EL FILE,GLOBAL=YES,TIMEPOINTS=T1")')
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,t,alp,bt,gmm,nb)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3),t(3,3),alp(4),bt(4),gmm(4)
!応力－ひずみ行列：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
!---固定軸回りのオイラー角？の回転行列
t=0.d0
t(1,1)=              cos(bt(nb))*cos(gmm(nb))
t(2,1)= sin(alp(nb))*sin(bt(nb))*cos(gmm(nb))+cos(alp(nb))*sin(gmm(nb))
t(3,1)=-cos(alp(nb))*sin(bt(nb))*cos(gmm(nb))+sin(alp(nb))*sin(gmm(nb))
t(1,2)=-cos(bt(nb))*sin(gmm(nb))                               
t(2,2)=-sin(alp(nb))*sin(bt(nb))*sin(gmm(nb))+cos(alp(nb))*cos(gmm(nb))
t(3,2)= cos(alp(nb))*sin(bt(nb))*sin(gmm(nb))+sin(alp(nb))*cos(gmm(nb))
t(1,3)= sin(bt(nb))                                                
t(2,3)=-sin(alp(nb))*cos(bt(nb))                                    
t(3,3)= cos(alp(nb))*cos(bt(nb))

return
end
