! 木材を直交異方性材料としてccxで解く場合の応力-ひずみ行列成分を
! 計算するプログラム。
implicit real*8(a-h,o-z)
dimension x(10000),y(10000),z(10000),es66(6,6),ew66(6,6)
! 木材の材料定数
ezz=6.d9 /1.d6
exx=ezz/25.d0
eyy=ezz/25.d0
! 上のヤング率と下のポアソン比は対称性を満たすように設定。
! exx=ezz/25だったら、poixy=poizx/25とか。
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数
gyz=ezz/15.d0
gxz=gyz
gxy=gyz
call zeromx(6,ew66)
 ew66(1,1)=  1.d0/exx ; ew66(1,2)=-poixy/exx ; ew66(1,3)=-poixz/exx
 ew66(2,1)=-poiyx/eyy ; ew66(2,2)=  1.d0/eyy ; ew66(2,3)=-poiyz/eyy
 ew66(3,1)=-poizx/ezz ; ew66(3,2)=-poizy/ezz ; ew66(3,3)=  1.d0/ezz
call gyaku(ew66)
 ew66(4,4)=gyz
 ew66(5,5)=gxz
 ew66(6,6)=gxy
!
!do i=1,6
! write(*,*) (ew66(i,j),j=1,6)
!end do
!
! 材料条件
!
 print'("*MATERIAL,NAME=mokuzai")'
 print'("*ELASTIC,TYPE=ORTHO")'
 print'(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a)',&
        ew66(1,1),',',ew66(1,2),',',ew66(2,2),',',ew66(1,3),',',&
       &ew66(2,3),',',ew66(3,3),',',ew66(6,6),',',ew66(5,5),','
 print*,ew66(4,4)
!
 print'("*SOLID SECTION,ELSET=Eall,MATERIAL=mokuzai")'
!
!
 end
!
!
! ひずみ−応力行列の逆数を取って、応力−ひずみ行列を求める。
!
subroutine gyaku(a)
 implicit real*8(a-h, o-z)
 dimension a(6,6), b(3,3)
  det=a(1,1)*a(2,2)*a(3,3)+a(2,1)*a(3,2)*a(1,3)&
& +a(3,1)*a(1,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)&
& -a(3,1)*a(2,2)*a(1,3)-a(1,1)*a(3,2)*a(2,3)
 b(1,1)=(-1.d0)**(1+1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))
 b(1,2)=(-1.d0)**(1+2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))
 b(1,3)=(-1.d0)**(1+3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2))
 b(2,1)=(-1.d0)**(2+1)*(a(1,2)*a(3,3)-a(3,2)*a(1,3))
 b(2,2)=(-1.d0)**(2+2)*(a(1,1)*a(3,3)-a(3,1)*a(1,3))
 b(2,3)=(-1.d0)**(2+3)*(a(1,1)*a(3,2)-a(3,1)*a(1,2))
 b(3,1)=(-1.d0)**(3+1)*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
 b(3,2)=(-1.d0)**(3+2)*(a(1,1)*a(2,3)-a(2,1)*a(1,3))
 b(3,3)=(-1.d0)**(3+3)*(a(1,1)*a(2,2)-a(2,1)*a(1,2))
do i=1,3
do j=1,3
 a(i,j)=b(j,i)/det
end do
end do
             return
             end
!
!
subroutine zeromx(n,a)
implicit real*8(a-h,o-z)
 dimension a(n,n)
  do i=1,n
  do j=1,n
     a(i,j)=0.d0
  end do
  end do
 return
end
!
!
