      program shyniten2 ver1 !2007/1/2
      implicit real*8(a-h,o-z)
      dimension nofn(8),b(3,3),e66(6,6),area(100,100),e77(6,6)
!
!######## 各断面の諸定数入力 ########
!      print(*,*)'x軸(断面の横の)の長さを入れてください'
!      read(*,*)zl
      zl=6.d-2
      zl=zl/2.d0
!
!      print*,'y軸(断面の縦の)の長さを入れてください'
!      read(*,*)zm
      zm=12.d-2
!
!       print*,'z軸(軸長)の長さを入れてください'
!       print*,'支点からはみ出している梁の長さも含めてください'
!       read*,ziku
!       ziku=ziku/2.d0
       ziku=1.18d0
!
!       print*,'支点からはみ出している長さを入れてください(片方だけ)'
!       read*,zikuhami
       zikuhami=1.d-1
!
!
!      print*,'上部に挿入する鋼板のx軸の長さを入れてください'
!      read(*,*)zlku
      zlku=9.d-3
      zlku=zlku/2.d0
!
!      print*,'上部に挿入する鋼板のy軸の長さを入れてください'
!      read(*,*)zmku
      zmku=4.d-2
!      print(*,*)'下部に挿入する鋼板のx軸の長さを入れて下さい'
!      read(*,*)zlkl
      zlkl=4.5d-3
!      
!      print*,'下部に挿入する鋼板のy軸の長さを入れて下さい'
!      read(*,*)zmkl
      zmkl=4.d-2
!
!      print*,'鋼板は図心軸からどのくらいのところに入れますか'
!      print*,'鋼板の厚みは考慮せず入れる位置だけ指定してください'
!      read(*,*) zlt
      zlt=1.d-2
      zlt=zlt-zlku-zlkl/2.d0
!
!######## 要素数入力 ########
!
!      print*,'x軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!      read(*,*) nx
      nx=20
      nx=nx/2
!
!      print*,'y軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!      read(*,*) ny
      ny=30
!
!       print*,'載荷する場所を入れて下さい'
!       print*,'ただし、梁中央からの距離を、梁の半分の軸長で割った値'
!       print*,'(支点からはみ出しているところを除く)
!       print*,'を通分して、分子,分母、1/3なら 1,3 と入れて下さい'
!       read*,ibunsi,ibunbo
       ibunsi=2
       ibunbo=9
!       print*,'今入力した分母を何倍かしてz方向の要素数にします'
!       print*,'何倍しますか(かけて120前後になるのが目安)'
!       print*,'支点からはみ出している部分の要素分割は'
!       print*,'梁の有効部分の要素分割にあわせて分割されます'
!       read*,ibunbokake
       ibunbokake=14
       nz0=ibunbokake*ibunbo
! 変数を整数系にした場合、右が実数系でも整数形になる
! また、少数は切り捨てられる 1.9=1 みたいに
       nzhami=zikuhami*nz0/ziku
       nz=nz0+nzhami
       nzd=ibunbokake*ibunsi
!
!
!      print*,'上部の鋼板だけのx軸方向の要素数(偶数)を入れてください'
!      read(*,*) nxku
      nxku=4
      nxku=nxku/2
!      print*,'上部の鋼板だけのy軸方向の要素数を入れてください'
!      read(*,*) nyku
      nyku=10
!
!      print*,'下部の鋼板だけのx軸方向の要素数を入れてください'
!      read(*,*) nxkl
      nxkl=1
!
!      print*,'下部の鋼板だけのy軸方向の要素数を入れてください'
!      read(*,*) nxkl
      nykl=10
!
!      print*,'鋼板の間の木材の要素数を入れて下さい'
!      read(*,*) nxt
      nxt=3
!
!
!######## 出力 #########
      open(7,file='shyniten2.inp')
!
!座標読み込み
      call zahyou(nx,ny,nz,nz0,nzhami,zl,zm,ziku,zikuhami,zlku,
     & zmku,zlkl,zmkl,nxku,nxkl,nyku,nykl,zlt,nnofn,zax,zay,zaz,nxt)
!
!要素番号読み込み
      call youso(nx,ny,nz,nxk,nyk,nofn,nofe,nxkl,nxku,nykl,nyku,
     &nxt)
!
!
!######## 境界条件 ########
!境界条件に必要な座標
!
!固定端の座標
       write(7,'("*NSET,NSET=FIX")')
       do i=1,(nx+1)*(ny+1)
       write(7,*)i,','
       end do
!図心を含むようにyz平面で切断した面の座標
      write(7,'("*NSET,NSET=SETUDAN")')
      do j=1,nz+1
      do i =1,ny+1 
      write(7,*)(j-1)*(ny+1)*(nx+1)+(i-1)*(nx+1)+1,','
      end do
      end do
!
!支点の面の一番下のx軸のy方向を拘束
      write(7,'("*NSET,NSET=SITEN")')
      do i=1,nx+1
      write(7,*)(nx+1)*(ny+1)*nz0+(nx+1)*ny+i,','
      end do
!
      write(7,'("*BOUNDARY")')
      write(7,'("FIX,3")')
      write(7,'("SETUDAN,1")')
      write(7,'("SITEN,2")')
!
!######## 材料の設定 ########
!******** 集成材の材料設定(圧縮) ********
      write(7,'("*MATERIAL,NAME=GLUP")')
      write(7,'("*ELASTIC,TYPE=ORTHOTROPIC")')
! 
!　ヤング率：exx,eyy,ezz
!      write(*,*)'z軸方向のヤング率を入れて下さい'
!      read(*,*) ezz
       ezz=9.553d3
       ezz=11.130d3
!      write(*,*)'x軸方向のヤング率を入れて下さい'
!      read(*,*) exx
       exx=ezz/25.d0
!      write(*,*)'y軸方向のヤング率を入れて下さい'
!      read(*,*) eyy
       eyy=ezz/25.d0
!
!　ポアソン比：poixy,poixz,poiyz
       poixy=0.016d0
       poixz=poixy
       poiyx=poixy
       poiyz=poixy
       poizx=0.4d0
       poizy=poizx
!
!　せん断弾性係数：gxy,gxz,gxy
       gxy=ezz/15.d0
       gxz=gxy
       gyz=gxy
!
      call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,
     & gxz,gxy,poiyx,poizx,poizy)
       write(7,'(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(7,'(f12.4,a,a)')e66(4,4),',','293'
!++++++++++ 塑性を入れる場合(集成材・圧縮) +++++++++
!鋼板挿入集成材梁の有限要素解析：
!千田知弘、後藤文彦、薄木征三、佐々木貴信：構造工学論文集Vol.53A  
!より集成材の降伏応力を引用
       write(7,'("*PLASTIC")')
!       write(7,'("26.366,0.")')
!       write(7,'("26.366,0.5")')
       write(7,'("32.55,0.")')
       write(7,'("32.55,0.5")')
!******** 集成材の材料設定(引張) ********
       write(7,'("*MATERIAL,NAME=GLUT")')
       write(7,'("*ELASTIC,TYPE=ORTHOTROPIC")')
       write(7,'(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(7,'(f12.4,a,a)')e66(4,4),',','293'
!******** 鋼材の材料設定 ********
      write(7,'("*MATERIAL,NAME=STE")')
      write(7,'("*ELASTIC,TYPE=ORTHO")')
      exxk=2.1d5
      eyyk=exxk
      ezzk=exxk
      poixyk=0.3d0
      poixzk=poixyk
      poiyxk=poixyk
      poiyzk=poixyk
      poizxk=poixyk
      poizyk=poixyk
      gxyk=8.076d4
      gxzk=gxyk
      gyzk=gxyk
      call ortho(e77,exxk,eyyk,ezzk,poixyk,poixzk,poiyzk,gyzk,
     & gxzk,gxyk,poiyxk,poizxk,poizyk)
       write(7,'(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a
     & ,f12.4,a)')
     &  e77(1,1),',',e77(1,2),',',e77(2,2),',',e77(1,3),',' 
     & ,e77(2,3),',',e77(3,3),',',e77(6,6),',',e77(5,5),','
       write(7,'(f12.4,a,a)')e77(4,4),',','293'
!+++++++++ 塑性を入れる場合(鋼材) +++++++
       write(7,'("*PLASTIC")')
       write(7,'("235.,0.")')
        write(7,'("235.,0.5")')
!
!#######################
! 要素に材料特性の付加 #
!#######################
! ここでどの要素がどんな材料か決定
      write(7,'("*SOLID SECTION,ELSET=glulamp,MATERIAL=GLUP")')
      write(7,'("*SOLID SECTION,ELSET=glulamt,MATERIAL=GLUT")')
      write(7,'("*SOLID SECTION,ELSET=ukou,MATERIAL=STE")')
      write(7,'("*SOLID SECTION,ELSET=lkou,MATERIAL=STE")')
!
! 
       write(7,'("*STEP")')
       write(7,'("*STATIC,SOLVER=SPOOLES")')
       write(7,'("*CLOAD")')
!
!載荷条件読み込み
      print*,'線載荷(x軸)にしますか、面載荷にしますか？'
      print*,'線載荷：0 , 面載荷：1'
      read*,isaikasentaku
      if(isaikasentaku.eq.0)then
      call line(nx,nxku,nxkl,nzd,nxt,zl,zlku,zlkl,zlt)
      else
      call menseki(nx,ny,nzd,nxk,zl,zlk,zm,zmk,nyk,nxku,nxkl,nyku,
     &nykl,nxt,zlku,zlkl,zmku,zmkl,nofnp,pakjyu,zakajyu,area,zlt)
      end if
!
       write(7,'("*NODE PRINT,NSET=Nall")')
       write(7,'("U")')
       write(7,'("*EL PRINT,ELSET=glulamp")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=glulamt")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=ukou")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=lkou")')
       write(7,'("S")')
       write(7,'("*NODE FILE")')
       write(7,'("U")')
       write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")')
       write(7,'("S")')
       write(7,'("*END STEP")')
!
      close(7)
      end
!###################### ここまでがメイン ###################
!
!######## 座標を出力するサブルーチン######
      subroutine zahyou(nx,ny,nz,nz0,nzhami,zl,zm,ziku,zikuhami,zlku,
     & zmku,zlkl,zmkl,nxku,nxkl,nyku,nykl,zlt,nnofn,zax,zay,zaz,nxt)
      implicit real*8(a-h,o-z)
       write(7,'("*NODE,NSET=Nall")')
       nnofn=0
       zykiten=zm/-2.d0
       do k=1,nz+1
       do j=1,ny+1
       do i=1,nx+1
! xの座標出力
        nnofn=nnofn+1
       if(i.le.2) then
        zax=(i-1)*zlku/(2*nxku-1)
       else if(i.le.nxku+1) then
        zax=zlku/(2*nxku-1)+(i-2)*zlku*2.d0/(2*nxku-1)
       else if(i.le.nxku+nxt+1) then
        zax=zlku+(i-nxku-1)*(zlt)/nxt
       else if(i.le.nxku+nxt+nxkl+1) then
        zax=zlku+zlt+(i-nxku-nxt-1)*zlku/nxkl
       else
        zax=zlku+zlt+zlkl+(i-nxku-nxt-nxkl-1)*(zl-zlku-zlt-zlkl)
     & /(nx-nxkl-nxku-nxt)
       end if
!
! yの座標出力
       if(j.le.nyku+1) then
        zay=zykiten+(j-1)*zmku/nyku
       else if(j.le.ny+1-nykl) then
        zay=zykiten+zmku+(j-nyku-1)*(zm-zmku-zmkl)/(ny-nyku-nykl)
       else
        zay=zykiten+zm-zmkl+(j-ny+nykl-1)*zmkl/nykl
       end if
!
! zの座標出力
       if(k.le.nz0+1)then
        zaz=(k-1)*(ziku-zikuhami)/nz0
       else
        zaz=ziku-zikuhami+(k-nz0-1)*zikuhami/nzhami
       end if
!
      write(7,'(i5,a,F10.7,a,F10.7,a,F10.7)')nnofn,',',zax,',',zay,
     & ',',zaz
       end do
       end do
       end do
!
       return
       end
!
!
!######## 要素を出力するサブルーチン ########
      subroutine youso(nx,ny,nz,nxk,nyk,nofn,nofe,nxkl,nxku,nykl,nyku,
     &nxt)
      implicit real*8(a-h,o-z)
      dimension nofn(8)
!cccccccccccc 鋼材と木材の節点番号出力 cccccccccccccccc
      do k=1,nz
      do j=1,ny
      do i=1,nx
!****** 上部鋼出力 ******
      if((i.le.nxku).and.(j.le.nyku))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=ukou")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
!****** 下部鋼出力 ******
       else if((i.ge.nxku+nxt+1).and.(i.le.nxku+nxt+nxkl).and.
     &(j.ge.ny-nykl+1))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=lkou")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!****** 集成材出力 ******
!++++++ 圧縮側 ++++++
       else if(j.le.ny/2)then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=glulamp")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
!++++++ 引張側 ++++++
       else 
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
      write(7,'("*ELEMENT,TYPE=C3D8,ELSET=glulamt")')
      write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
       end if
      end do
      end do
      end do
      return
      end
!
!
!######## 各要素の面積に応じて載荷するサブルーチン ########
      subroutine menseki(nx,ny,nzd,nxk,zl,zlk,zm,zmk,nyk,nxku,nxkl,nyku,
     &nykl,nxt,zlku,zlkl,zmku,zmkl,nofnp,pakjyu,zakajyu,area,zlt)
      implicit real*8(a-h,o-z)
      dimension area(100,100),yoko(100),tate(100),ax(100),ay(100)
!******* 各要素の面積を計算 ******
       do j=1,ny+1
       do i=1,nx+1
! xの座標出力
       if(i.le.2) then
        yoko(i)=(i-1)*zlku/(2*nxku-1)
       else if(i.le.nxku+1) then
        yoko(i)=zlku/(2*nxku-1)+(i-2)*zlku*2.d0/(2*nxku-1)
       else if(i.le.nxku+nxt+1) then
        yoko(i)=zlku+(i-nxku-1)*(zlt)/nxt
       else if(i.le.nxku+nxt+nxkl+1) then
        yoko(i)=zlku+zlt+(i-nxku-nxt-1)*zlku/nxkl
       else
        yoko(i)=zlku+zlt+zlkl+(i-nxku-nxt-nxkl-1)*(zl-zlku-zlt-zlkl)
     & /(nx-nxkl-nxku-nxt)
       end if
!
! yの座標出力
       if(j.le.nyku+1) then
        tate(j)=(j-1)*zmku/nyku
       else if(j.le.ny+1-nykl) then
        tate(j)=zmku+(j-nyku-1)*(zm-zmku-zmkl)/(ny-nyku-nykl)
       else
        tate(j)=zm-zmkl+(j-ny+nykl-1)*zmkl/nykl
       end if
!
       end do
       end do
!
! 各要素の横の長さ
      yoko(1)=0
      do i=1,nx
      ax(i)=yoko(i+1)-yoko(i)
      end do
! 各要素の縦の長さ
      tate(1)=0
      do j=1,ny
      ay(j)=tate(j+1)-tate(j)
      end do
! 各要素の面積
      do j=1,ny
      do i=1,nx
      area(j,i)=ay(j)*ax(i)
      end do
      end do
!****** 自由端の各節点に載荷 ******
      write(*,'("載荷するトータルの荷重を入力して下さい")')
      read*,zkajyu
      zakajyu=zkajyu/zl/zm
! 変位を正で出したいとき、つまり、
! 座標のプラス方向に載荷する場合は正、その逆には負(-1倍)
!      zakajyu=zakajyu*-1.d0
!片持ち梁だから2で割る、さらに半解析だから2で割る
      zakajyu=zakajyu/2.d0/2.d0
!++++++ 4要素と接する節点 ++++++
! 対称面を含む要素の節点
      do j=2,ny
       pakajyu=(2.d0*area(j-1,1)+area(j-1,2)+2.d0*area(j,1)+area(j,2))
     &/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nzd+(j-1)*(nx+1)+2
       write(7,*)nofnp,',','2',',',pakajyu
      end do
! その他 
      do j=2,ny
      do i=3,nx
       pakajyu=(area(j-1,i-1)+area(j-1,i)+area(j,i-1)+area(j,i))/4.d0
     &*zakajyu
       if((i.gt.1).and.(i.lt.nx+1).and.(j.gt.1).and.(j.lt.ny+1))then
       nofnp=(nx+1)*(ny+1)*nzd+(j-1)*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
      end if
      end do
      end do
!
!++++++ 2要素と接する節点 ++++++
! 上(対称面を含む要素)
       pakajyu=(2.d0*area(1,1)+area(1,2))/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nzd+2
       write(7,*)nofnp,',','2',',',pakajyu
! 上
      do i =3,nx
       pakajyu=(area(1,i-1)+area(1,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nzd+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 下(対称面を含む要素)
       pakajyu=(2.d0*area(ny,1)+area(ny,2))/4.d0*zakajyu
       nofnp=(nx+1)*(ny+1)*nzd+ny*(nx+1)+2
       write(7,*)nofnp,',','2',',',pakajyu
! 下
      do i =3,nx
       pakajyu=(area(ny,i-1)+area(ny,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nzd+ny*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 右
      do j =2,ny
       pakajyu=(area(j-1,nx)+area(j,nx))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nzd+j*(nx+1)
      write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
!++++++ どことも接しない節点 ++++++
!
      write(7,*)(nx+1)*(ny+1)*nzd+nx+1,',','2',',',
     &area(1,nx)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nzd+ny*(nx+1)+nx+1,',','2',',',
     &area(ny,nx)/4.d0*zakajyu
!
      return
      end
! 
!
!######## 各要素のx方向の長さに応じて載荷するサブルーチン ########
!
      subroutine line(nx,nxku,nxkl,nzd,nxt,zl,zlku,zlkl,zlt)
      implicit real*8(a-h,o-z)
      dimension area(100,100),yoko(100),tate(100),ax(100),ay(100)
!******* 各要素のx軸方向の長さを計算 ******
       do i=1,nx+1
! xの座標出力
       if(i.le.2) then
        yoko(i)=(i-1)*zlku/(2*nxku-1)
       else if(i.le.nxku+1) then
        yoko(i)=zlku/(2*nxku-1)+(i-2)*zlku*2.d0/(2*nxku-1)
       else if(i.le.nxku+nxt+1) then
        yoko(i)=zlku+(i-nxku-1)*(zlt)/nxt
       else if(i.le.nxku+nxt+nxkl+1) then
        yoko(i)=zlku+zlt+(i-nxku-nxt-1)*zlku/nxkl
       else
        yoko(i)=zlku+zlt+zlkl+(i-nxku-nxt-nxkl-1)*(zl-zlku-zlt-zlkl)
     & /(nx-nxkl-nxku-nxt)
       end if
      end do
! 各要素の横の長さ
      yoko(1)=0
      do i=1,nx
      ax(i)=yoko(i+1)-yoko(i)
      end do
!
!****** 各節点に載荷 ******
      write(*,'("載荷するトータルの荷重を入力して下さい")')
      read*,zkajyu
      zakajyu=zkajyu/zl
! 変位を正で出したいとき、つまり、
! 座標のプラス方向に載荷する場合は正、その逆には負(-1倍)
!      zakajyu=zakajyu*-1.d0
!片持ち梁だから2で割る、さらに半解析だから2で割る
      zakajyu=zakajyu/2.d0/2.d0
!
! 対称面のすぐ横の節点
      write(7,*)(nx+1)*(ny+1)*nzd+2,',','2',',',(ax(1)+ax(2)/2.d0)
     &*zakajyu
!
! 対称面から2節点と反対側の1節点を除いた節点
      do i=3,nx
      pakajyu=(ax(i-1)+ax(i))/2.d0*zakajyu
      nofnp=(nx+1)*(ny+1)*nzd+i
      write(7,*)nofnp,',','2',',',pakajyu
      end do
!
! 対称面から一番離れたはじっこの節点
      write(7,*)(nx+1)*(ny+1)*nzd+nx+1,',','2',',',ax(nx)*zakajyu/2.d0
!
      return
      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
!
!****** せん断補正係数を入力 ******
!      write(*,*)'xy面のせん断弾性係数を入れて下さい'
!      read(*,*) gxy
      e66(4,4)=gyz
!      write(*,*)'xz面のせん断弾性係数を入れて下さい'
!      read(*,*) gxz
      e66(5,5)=gxz
!      write(*,*)'yz面のせん断弾性係数を入れて下さい'
!      read(*,*) gyz
      e66(6,6)=gxy
!
!
      return
      end

