! まずは、beam8p.inpと同じデータを、要素番号と節点番号の
! 並び順を変えて任意の要素数で生成するプログラムを作ってみる。
implicit real*8(a-h,o-z)
dimension x(10000),y(10000),z(10000),es66(6,6),ew66(6,6)
! 節点座標：x(節点番号),y(節点番号),z(節点番号)
 print'("**")'
 print'("** beam8p.inpと同じデータを要素番号と節点番号の")'
 print'("** 並び順(と要素数)を変えたデータ")'
 print'("**")'
 print'("*NODE, NSET=Nall")'
! 桁幅(x方向)：b(m)
bs=6.d-3         ! 鋼板の桁幅(x方向)：bs(m)
bw=840.d-3/2     ! 木材の桁幅(x方向)：bw(m)
hw=120.d-3       ! 木材の桁高(y方向)：hw(m)
h=500.d-3        ! 桁高(y方向)：h(m)
hs=h-hw*2        ! 鋼板中央部の高さ(m)
ell=3.d0         ! 梁の軸長(z方向)：ell(m)
p=0.1d0 !中央集中荷重(MN)
!
!bs=2.d-3: bw=8.d-3 hw=4.d-3: h=10.d-3: hs=3.d-3: ell=10.d-3
!
! x,y,z方向の要素分割数：nx,ny,nz
nx=2
ny=4
nz=4
nxs=1  !鋼板の幅(x)方向の要素分割数
nxw=nx-nxs !木材の幅(x)方向の要素分割数
nyw=1  !木材の高さ(y)方向の要素分割数
nys=ny-nyw*2  !鋼板中央部の高さ(y)方向の要素分割数
nsiten=1 !梁両端から支点までの要素数
nmen=(nx+1)*(nyw+1)*2+(nxs+1)*(nys-1) !断面の全節点数
!
! 1要素のx,y,z方向の長さ：xn,yn,zn
xw=bw/nxw !木材要素1要素のx方向の長さ
xs=bs/nxs !木材要素1要素のx方向の長さ
yw=hw/nyw !木材要素1要素のy方向の長さ
ys=hs/nys !鋼板中央部要素1要素のy方向の長さ
zn=ell/nz !各要素の1要素のz方向の長さ
!
! print*,'要素数：',(2*nx*nyw+nxs*nys)*nz
! print*,'節点数：',(2*(nx+1)*(nyw+1)+(nxs+1)*(nys-1))*(nz+1)
!
! 鋼板の材料定数
es=206.d9 /1.d6
pois=0.3d0
gs=es/2.d0/(1.d0+pois)
call zeromx(6,ew66)
 es66(1,1)= 1.d0/es ; es66(1,2)=-pois/es ; es66(1,3)=-pois/es
 es66(2,1)=-pois/es ; es66(2,2)= 1.d0/es ; es66(2,3)=-pois/es
 es66(3,1)=-pois/es ; es66(3,2)=-pois/es ; es66(3,3)= 1.d0/es
 es66(4,4)=gs !逆行列を取るのは左上の3*3だけなので
 es66(5,5)=gs !右下3*3は最初からGを入れておく
 es66(6,6)=gs
call gyaku(es66)
! 木材の材料定数
ezz=6.d9 /1.d6
exx=ezz/25.d0
eyy=ezz/25.d0
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
!
!
! 節点番号ごとの節点座標
! 節点番号：nset
   nset=0
do k=0,nz
do j=0,ny
 if(j<=nyw) then
  do i=0,nx
   nset=nset+1
   if(i<=nxs) then 
    x(nset)=xs*i
    y(nset)=yw*j
    z(nset)=zn*k
    print*,nset,',',x(nset),',',y(nset),',',z(nset)
   else
    x(nset)=xs*nxs+xw*(i-nxs)
    y(nset)=yw*j
    z(nset)=zn*k
    print*,nset,',',x(nset),',',y(nset),',',z(nset)
   end if
  end do
 else if(j<ny-nyw) then
  do i=0,nxs
   nset=nset+1
    x(nset)=xs*i
    y(nset)=yw*nyw+ys*(j-nyw)
    z(nset)=zn*k
    print*,nset,',',x(nset),',',y(nset),',',z(nset)
  end do
 else
  do i=0,nx
   nset=nset+1
   if(i<=nxs) then 
    x(nset)=xs*i
    y(nset)=yw*nyw+ys*nys+yw*(j-nyw-nys)
    z(nset)=zn*k
    print*,nset,',',x(nset),',',y(nset),',',z(nset)
   else
    x(nset)=xs*nxs+xw*(i-nxs)
    y(nset)=yw*nyw+ys*nys+yw*(j-nyw-nys)
    z(nset)=zn*k
    print*,nset,',',x(nset),',',y(nset),',',z(nset)
   end if
  end do
 end if
end do
end do
!
!
! 要素番号とその要素の8節点の節点番号
 nyou=0
do k=1,nz
do j=1,ny
 if(j<=nyw) then
  do i=1,nx
   nyou=nyou+1
   n1=nmen*(k-1)+(nx+1)*(j-1)+i
   n2=n1+1
   n3=n1+(nx+1)
   n4=n3+1
   n5=n1+nmen
   n6=n2+nmen
   n7=n3+nmen
   n8=n4+nmen
   if(i<=nxs) then
    print'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")'
   else
    print'("*ELEMENT, TYPE=C3D8, ELSET=mokuzai")'
   end if
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
     &   nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
  end do
 else if(j==nyw+1) then
  do i=1,nxs
   nyou=nyou+1
   n1=nmen*(k-1)+(nx+1)*nyw+i
   n2=n1+1
   n3=n1+(nx+1)
   n4=n3+1
   n5=n1+nmen
   n6=n2+nmen
   n7=n3+nmen
   n8=n4+nmen
   if(i<=nxs) then
    print'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")'
   else
    print'("*ELEMENT, TYPE=C3D8, ELSET=mokuzai")'
   end if
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
     &   nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
  end do
 else if(j<=ny-nyw) then
  do i=1,nxs
   nyou=nyou+1
   n1=nmen*(k-1)+(nx+1)*(nyw+1)+(nxs+1)*(j-nyw-2)+i
   n2=n1+1
   n3=n1+(nxs+1)
   n4=n3+1
   n5=n1+nmen
   n6=n2+nmen
   n7=n3+nmen
   n8=n4+nmen
   if(i<=nxs) then
    print'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")'
   else
    print'("*ELEMENT, TYPE=C3D8, ELSET=mokuzai")'
   end if
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
     &   nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
  end do
 else !ny-nyw<j<ny
  do i=1,nx
   nyou=nyou+1
   n1=nmen*(k-1)+(nx+1)*(nyw+1)+(nxs+1)*(nys-1)+(nx+1)*(j-nyw-nys-1)+i
   n2=n1+1
   n3=n1+(nx+1)
   n4=n3+1
   n5=n1+nmen
   n6=n2+nmen
   n7=n3+nmen
   n8=n4+nmen
   if(i<=nxs) then
    print'("*ELEMENT, TYPE=C3D8, ELSET=kouhan")'
   else
    print'("*ELEMENT, TYPE=C3D8, ELSET=mokuzai")'
   end if
   print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
     &   nyou,',',n1,',',n3,',',n7,',',n5,',',n2,',',n4,',',n8,',',n6
  end do
 end if
end do
end do
!
!
! 境界条件
! 節点ごとに異なる拘束条件を与える場合をふまえて、
! 拘束条件の種類に対して適当な名前(FIXとか)を付ける
!
! 幅員中央の対称面を拘束
print'("*NSET,NSET=taisyou")'
do k=1,nz+1
do j=1,ny+1
 if(j<=nyw+1) then
  print'(i6,a)',nmen*(k-1)+(nx+1)*j,','
 else if (j>ny-nyw) then
  print'(i6,a)',nmen*(k-1)+(nx+1)*(nyw+1)+(nxs+1)*(nys-1)+(nx+1)*(j-nyw-nys),','
 end if
end do
end do
!
! ヒンジ支承
print'("*NSET,NSET=hinzi")'
nowari=nmen*(nsiten+1)
nhazime=nowari-nx
do i=nhazime,nowari
 print'(i6,a)',i,','
end do
! ローラー支承
print'("*NSET,NSET=rooraa")'
nowari=nset-nmen*nsiten
nhazime=nowari-nx
do i=nhazime,nowari
 print'(i6,a)',i,','
end do
!
! 上記(NSET=taisyou)の拘束節点の拘束変位
 print'("*BOUNDARY")'
 print'("taisyou,1")'
 print'("hinzi,2,3")'
 print'("rooraa,2")'
!
! 材料条件
!
 print'("*MATERIAL,NAME=kouhan")'
 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)',&
        es66(1,1),',',es66(1,2),',',es66(2,2),',',es66(1,3),',',&
       &es66(2,3),',',es66(3,3),',',es66(6,6),',',es66(5,5),','
 print*,es66(4,4),',',293.d0
!
!
 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),',',293.d0
!
 print'("*SOLID SECTION,ELSET=mokuzai,MATERIAL=mokuzai")'
 print'("*SOLID SECTION,ELSET=kouhan,MATERIAL=kouhan")'
!
! 載荷条件
print'("*NSET,NSET=saika")'
nhazime=nmen*(nz/2)+1
nowari=nhazime+nx
do i=nhazime,nowari
 print'(i6,a)',i,','
end do
!
 print'("*STEP")'
 print'("*STATIC,SOLVER=SPOOLES")'
 print'("*CLOAD")'
!  y方向に載荷
 pn=p/(nx+1)
 write(*,*)'saika,2,',pn
!
 print'("*NODE PRINT,NSET=Nall")'
 print'("U")'
 print'("*EL PRINT,ELSET=kouhan")'
 print'("S")'
 print'("*EL PRINT,ELSET=mokuzai")'
 print'("S")'
 print'("*NODE FILE")'
 print'("U")'
 print'("*EL FILE,POSITION=AVERAGED AT NODES")'
 print'("S")'
 print'("*END STEP")'
!
!
!
 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
!
!
