implicit real*8(a-h, o-z)
dimension f(30),& !節点力ベクトル
& d(30),& !節点変位ベクトル
& x(30),& !境界条件（拘束節点に0, その他に1が入る)
& ss(30,30),& !全体剛性行列
& xv(15),xw(15),& !節点ごとの変位境界条件
& fy(15),fz(15),& !節点ごとの荷重条件
& idari(15),migi(15),& !各要素の左節点番号、右節点番号
& s(15,4,4),& !要素剛性行列
& th(15),& !各要素の回転角
& t(4,4),& !座標変換行列（1要素ごとに計算）
& s44(4,4),& !要素剛性行列（1要素ごとに移し替える入れ物）
& ts(4,4),& !TK（1要素ごとに計算）
& tst(4,4) !TKT（1要素ごとに計算）
dimension sn(4) !節点力の出力のために追加

pi=asin(1.d0)*2.d0
nyou=15 !要素数
nset=9 !節点数

! 初期化
do n=1,nyou
do i=1,4
do j=1,4
s(n,i,j)=0.d0
end do
end do
end do

do i=1,30
d(i)=0.d0 !; f(i)=0.d0; x(i)=0.d0
do j=1,30
ss(i,j)=0.d0
end do
end do

do i=1,15
xv(i)=1.d0; xw(i)=1.d0 !境界条件は1で初期化
fy(i)=0.d0; fz(i)=0.d0
end do

E=14.5d3!MPa
ell=250.d0!mm
a=12.d0**2.d0!mm^2
sk1=E*a/ell
! 要素1の剛性マトリクス
s(1,2,2)= sk1; s(1,2,4)=-sk1
s(1,4,2)=-sk1; s(1,4,4)= sk1
! 要素2~15の剛性マトリクス
sk2=sk1
do k=2,15
do j=1,4
do i=1,4
s(k,j,i)=s(1,j,i)
enddo
enddo
enddo

!各要素の左節点番号、右節点番号
idari(1)=1; migi(1)=2
idari(2)=2; migi(2)=3
idari(3)=3; migi(3)=4
idari(4)=4; migi(4)=5
idari(5)=1; migi(5)=6
idari(6)=2; migi(6)=6
idari(7)=2; migi(7)=7
idari(8)=3; migi(8)=7
idari(9)=3; migi(9)=8
idari(10)=4; migi(10)=8
idari(11)=4; migi(11)=9
idari(12)=5; migi(12)=9
idari(13)=6; migi(13)=7
idari(14)=7; migi(14)=8
idari(15)=8; migi(15)=9

!各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th(1)=0.d0
th(2)=0.d0
th(3)=0.d0
th(4)=0.d0
th(5)=pi/3
th(6)=pi*2.d0/3
th(7)=pi/3
th(8)=pi*2.d0/3
th(9)=pi/3
th(10)=pi*2.d0/3
th(11)=pi/3
th(12)=pi*2.d0/3
th(13)=0.d0
th(14)=0.d0
th(15)=0.d0

!境界条件
xv(1)=0.d0; xw(1)=0.d0
xv(5)=0.d0!; xw(3)=0.d0

do i=1,15
x(2*i-1)=xv(i)
x(2*i)=xw(i)
end do
!
!載荷条件
fy(3)=300!N
!
do i=1,15
f(2*i-1)=fy(i)
f(2*i)=fz(i)
end do
!
!
!要素ごとのTKTの計算
do n=1,nyou
!
! 座標変換マトリクス
!
call zahyou(t,th(n))
!
do i=1,4
do j=1,4
s44(i,j)=s(n,i,j)
end do
end do
!
call mxtmx(t,s44,ts)
!
call mxmx(ts,t,tst)
!!
!!重ねあわせ
i1=2*idari(n)-1
i2=2*idari(n)
m1=2*migi(n)-1
m2=2*migi(n)
!
ss(i1,i1)=ss(i1,i1)+tst(1,1)
ss(i1,i2)=ss(i1,i2)+tst(1,2)
ss(i1,m1)=ss(i1,m1)+tst(1,3)
ss(i1,m2)=ss(i1,m2)+tst(1,4)
!         !
ss(i2,i1)=ss(i2,i1)+tst(2,1)
ss(i2,i2)=ss(i2,i2)+tst(2,2)
ss(i2,m1)=ss(i2,m1)+tst(2,3)
ss(i2,m2)=ss(i2,m2)+tst(2,4)
!         !
ss(m1,i1)=ss(m1,i1)+tst(3,1)
ss(m1,i2)=ss(m1,i2)+tst(3,2)
ss(m1,m1)=ss(m1,m1)+tst(3,3)
ss(m1,m2)=ss(m1,m2)+tst(3,4)
!         !
ss(m2,i1)=ss(m2,i1)+tst(4,1)
ss(m2,i2)=ss(m2,i2)+tst(4,2)
ss(m2,m1)=ss(m2,m1)+tst(4,3)
ss(m2,m2)=ss(m2,m2)+tst(4,4)
!
!
end do !n要素について
!!
!!
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!print*
!!
!　境界条件を入れる
do i=1,nset*2
do j=1,nset*2
  ss(i,j)=x(i)*ss(i,j)
  ss(j,i)=x(i)*ss(j,i)
end do
end do
do i=1,nset*2
if(x(i)<1.d-3) then
    ss(i,i)=1.d0
end if
end do
!!
!!
!
!print*,'f=',(f(j),j=1,18)
!
!
call gausu(nset*2,ss,d,f)
!
!print*
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!
do n=1,nset
print*,'節点番号：',n,'v=',d(n*2-1),'w=',d(n*2)
end do
!追加
   !要素ごとのTKTの計算
   do n=1,nyou
   !
   ! 座標変換マトリクス
   !
   call zahyou(t,th(n))
   !
   do i=1,4
   do j=1,4
   s44(i,j)=s(n,i,j)
   end do
   end do
   !
   call mxtmx(t,s44,ts)
   !
   call mxmx(ts,t,tst)
   !!
    sn=0.0
    do i=1,4
      do j=1,2
         sn(i)=sn(i)+tst(i,j)*d((idari(n)-1)*2+j)
      enddo
      do j=3,4
         sn(i)=sn(i)+tst(i,j)*d((migi(n)-1)*2+j-2)
      enddo
    enddo
    print'(A,I2,A,I2,A,f15.10,A,I2,A,f15.10)','部材番号： ',n,&
        &'  S(',idari(n),')=',sn(1),'  ,N(',idari(n),')=',sn(2)
    print'(A,A,I2,A,f15.10,A,I2,A,f15.10)','             ',&
        &'  S(',migi(n),')=',sn(3),'  ,N(',migi(n),')=',sn(4)
   !
   end do !n要素について
!ここまで追加
!print*,'d=',(d(j),j=1,18)
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
! 以上がメインプログラム
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
!
subroutine mxmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(4,4), b(4,4), c(4,4)
!
!
do i=1,4
do j=1,4
c(i,j)=0.
 do k=1,4
  c(i,j)=c(i,j)+a(i,k)*b(k,j)
 end do
end do
end do
!
return
end
!
subroutine zahyou(t,theta)
implicit real*8(a-h, o-z)
dimension t(4,4)
do i=1,2
do j=1,2
 t(i+2,j)=0.d0
 t(i,j+2)=0.d0
end do
end do
t(1,1)= cos(theta)
t(1,2)= sin(theta)
t(2,1)=-sin(theta)
t(2,2)= cos(theta)
t(3,3)=t(1,1)
t(3,4)=t(1,2)
t(4,3)=t(2,1)
t(4,4)=t(2,2)
return
end
!
subroutine mxtmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(4,4), b(4,4), c(4,4)
!
do i=1,4
do j=1,4
c(i,j)=0.
 do k=1,4
  c(i,j)=c(i,j)+a(k,i)*b(k,j)
 end do
end do
end do
!
!
return
end
!
!
subroutine gausu(n,a,x,b)
implicit real*8(a-h, o-z)
!dimension a(n,n),x(n),b(n)
dimension a(30,30),x(30),b(30)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15
!!
!print*
!do i=1,8
!print'(8f10.3)', (a(i,j),j=1,8)
!end do
!
!
do k=1,n-1  !a(k,k)を消去
do i=k+1,n  !k+1行からn行まで
do j=k+1,n    !k+1列からn列まで
   a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k)
end do
   b(i)=b(i)-b(k)*a(i,k)/a(k,k)
end do
end do
!
! 後退代入
x(n)=b(n)/a(n,n)
do k=n,1,-1
   akjxj=0.d0
do j=k+1,n
   akjxj=akjxj+a(k,j)*x(j)
end do
   x(k)=(b(k)-akjxj)/a(k,k)
end do
!
return
end
!
!

