implicit real*8(a-h, o-z)
dimension f(45),& !節点力ベクトル
& d(45),& !節点変位ベクトル
& x(45),& !境界条件（拘束節点に0, その他に1が入る)
& ss(45,45),& !全体剛性行列
& xv(15),xw(15),xt(15),& !節点ごとの変位境界条件:v,w,θ
& fy(15),fz(15),cx(15),& !節点ごとの荷重条件:S,N,M
& idari(15),migi(15),& !各要素の左節点番号、右節点番号
& s(15,6,6),& !要素剛性行列
& th(15),& !各要素の回転角
& t(6,6),& !座標変換行列（1要素ごとに計算）
& s66(6,6),& !要素剛性行列（1要素ごとに移し替える入れ物）
& ts(6,6),& !TK（1要素ごとに計算）
& tst(6,6) !TKT（1要素ごとに計算）
dimension snm(6) !節点力の出力のために追加
!
pi=asin(1.d0)*2.d0
!
nyou=15 !要素数
nset=9 !節点数
!
! 初期化
do n=1,nyou
do i=1,6
do j=1,6
s(n,i,j)=0.d0
end do
end do
end do
!
do i=1,45
d(i)=0.d0
do j=1,45
ss(i,j)=0.d0
end do
end do
!
do i=1,15
xv(i)=1.d0; xw(i)=1.d0; xt(i)=1.d0 !境界条件は1で初期化
fy(i)=0.d0; fz(i)=0.d0; cx(i)=0.d0
end do
!
! 例題1
! 要素1の剛性マトリクス
el=250.d0 !1本のL[mm]
area=12.d0**2.d0 !A[mm^2]
dan2=12.d0*12.d0**3/12 !I[mm^4]
young=14.5d3 !E[N/mm^2]
poi=0.4 !ポアソン比
hosei=10.d0*(1.d0+poi)/(12.d0+11.d0*poi) !矩形断面のせん断補正係数
ea=young*area !EA[N]
ei=young*dan2 !EI(N/mm^2・mm^4)
g=young/15 !木材のせん断弾性係数
rmd=el/sqrt(dan2/(hosei*area))
alp=young/(g*(rmd)**2)
print*,"g=",g
print*,"k=",hosei
print*,"rmd=",rmd
print*,"alp=",alp
s(1,2,2)=       ea/el                    ;s(1,2,5)=                  -ea/el
s(1,5,2)=      -ea/el                    ;s(1,5,5)=                   ea/el
s(1,1,1)= 12.d0*ei/el**3/(1.d0+12.d0*alp);s(1,1,3)=             -6.d0*ei/el**2/(1.d0+12.d0*alp)
s(1,1,4)=-12.d0*ei/el**3/(1.d0+12.d0*alp);s(1,1,6)=             -6.d0*ei/el**2/(1.d0+12.d0*alp)
s(1,3,1)= -6.d0*ei/el**2/(1.d0+12.d0*alp);s(1,3,3)=  (4.d0+12.d0*alp)*ei/el/(1.d0+12.d0*alp)
s(1,3,4)=  6.d0*ei/el**2/(1.d0+12.d0*alp);s(1,3,6)=  (2.d0-12.d0*alp)*ei/el/(1.d0+12.d0*alp)
s(1,4,1)=-12.d0*ei/el**3/(1.d0+12.d0*alp);s(1,4,3)=              6.d0*ei/el**2/(1.d0+12.d0*alp)
s(1,4,4)= 12.d0*ei/el**3/(1.d0+12.d0*alp);s(1,4,6)=              6.d0*ei/el**2/(1.d0+12.d0*alp)
s(1,6,1)= -6.d0*ei/el**2/(1.d0+12.d0*alp);s(1,6,3)=  (2.d0-12.d0*alp)*ei/el/(1.d0+12.d0*alp)
s(1,6,4)=  6.d0*ei/el**2/(1.d0+12.d0*alp);s(1,6,6)=  (4.d0+12.d0*alp)*ei/el/(1.d0+12.d0*alp)
! 要素2~15の剛性マトリクス
do k=2,15
do i=1,6
do j=1,6
   s(k,i,j)=s(1,i,j)
end do
end do
end do
!
!do i=1,6
!print'(9f10.2)',(s(1,i,j),j=1,9)
!end do
!print*
!!
!
!
!各要素の左節点番号、右節点番号
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!; xt(4)=0.d0
xv(5)=0.d0!; xw(3)=0.d0!; xt(3)=0.d0
!
!
do i=1,15
x(3*i-2)=xv(i)
x(3*i-1)=xw(i)
x(3*i)=xt(i)
end do
!
!載荷条件
!cx(2)=1.d0 ! 例題1
fy(3)=300!N
!
do i=1,15
f(3*i-2)=fy(i)
f(3*i-1)=fz(i)
f(3*i)=cx(i)
end do
!
!
!要素ごとのTKTの計算
do n=1,nyou
!
! 座標変換マトリクス
!
call zahyou(t,th(n))
!
do i=1,6
do j=1,6
s66(i,j)=s(n,i,j)
end do
end do
!
call mxtmx(t,s66,ts)
!
call mxmx(ts,t,tst)
!
!do i=1,6
!print'(9f10.2)',(tst(i,j),j=1,9)
!end do
!print*
!!!
!!重ねあわせ
i1=3*idari(n)-2
i2=3*idari(n)-1
i3=3*idari(n)
m1=3*migi(n)-2
m2=3*migi(n)-1
m3=3*migi(n)
!
ss(i1,i1)=ss(i1,i1)+tst(1,1)
ss(i1,i2)=ss(i1,i2)+tst(1,2)
ss(i1,i3)=ss(i1,i3)+tst(1,3)
ss(i1,m1)=ss(i1,m1)+tst(1,4)
ss(i1,m2)=ss(i1,m2)+tst(1,5)
ss(i1,m3)=ss(i1,m3)+tst(1,6)
!
ss(i2,i1)=ss(i2,i1)+tst(2,1)
ss(i2,i2)=ss(i2,i2)+tst(2,2)
ss(i2,i3)=ss(i2,i3)+tst(2,3)
ss(i2,m1)=ss(i2,m1)+tst(2,4)
ss(i2,m2)=ss(i2,m2)+tst(2,5)
ss(i2,m3)=ss(i2,m3)+tst(2,6)
!
ss(i3,i1)=ss(i3,i1)+tst(3,1)
ss(i3,i2)=ss(i3,i2)+tst(3,2)
ss(i3,i3)=ss(i3,i3)+tst(3,3)
ss(i3,m1)=ss(i3,m1)+tst(3,4)
ss(i3,m2)=ss(i3,m2)+tst(3,5)
ss(i3,m3)=ss(i3,m3)+tst(3,6)
!         !
!
ss(m1,i1)=ss(m1,i1)+tst(4,1)
ss(m1,i2)=ss(m1,i2)+tst(4,2)
ss(m1,i3)=ss(m1,i3)+tst(4,3)
ss(m1,m1)=ss(m1,m1)+tst(4,4)
ss(m1,m2)=ss(m1,m2)+tst(4,5)
ss(m1,m3)=ss(m1,m3)+tst(4,6)
!
ss(m2,i1)=ss(m2,i1)+tst(5,1)
ss(m2,i2)=ss(m2,i2)+tst(5,2)
ss(m2,i3)=ss(m2,i3)+tst(5,3)
ss(m2,m1)=ss(m2,m1)+tst(5,4)
ss(m2,m2)=ss(m2,m2)+tst(5,5)
ss(m2,m3)=ss(m2,m3)+tst(5,6)
!
ss(m3,i1)=ss(m3,i1)+tst(6,1)
ss(m3,i2)=ss(m3,i2)+tst(6,2)
ss(m3,i3)=ss(m3,i3)+tst(6,3)
ss(m3,m1)=ss(m3,m1)+tst(6,4)
ss(m3,m2)=ss(m3,m2)+tst(6,5)
ss(m3,m3)=ss(m3,m3)+tst(6,6)
!
!
end do !n要素について
!!
!!
!do i=1,8
!print'(8f10.3)', (ss(i,j),j=1,8)
!end do
!print*
!!
!　境界条件を入れる
do i=1,nset*3
do j=1,nset*3
  ss(i,j)=x(i)*ss(i,j)
  ss(j,i)=x(i)*ss(j,i)
end do
end do
do i=1,nset*3
if(x(i)<1.d-3) then
    ss(i,i)=1.d0
end if
end do
!!
!!
!do i=1,9
!print'(9f10.2)',(ss(i,j),j=1,9)
!end do
!
!print*,'f=',(f(j),j=1,18)
!
!
call gausu(nset*3,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*3-2),'w=',d(n*3-1), 'th=',d(n*3)
end do
!追加
    print'(A)','                                                    曲げモーメントの単位に注意!'
   !要素ごとのTKTの計算
   do n=1,nyou
   !
   ! 座標変換マトリクス
   !
   call zahyou(t,th(n))
   !
   do i=1,6
   do j=1,6
   s66(i,j)=s(n,i,j)
   end do
   end do
   !
   call mxtmx(t,s66,ts)
   !
   call mxmx(ts,t,tst)
   !
    snm=0.0
    do i=1,6
      do j=1,3
         snm(i)=snm(i)+tst(i,j)*d((idari(n)-1)*3+j)
      enddo
      do j=4,6
         snm(i)=snm(i)+tst(i,j)*d((migi(n)-1)*3+j-3)
      enddo
    enddo
    print'(A,I2,A,I2,A,f15.10,A,I2,A,f15.10,A,I2,A,f15.10)','部材番号： ',n,&
        &'   S(',idari(n),')=',snm(1),'  ,N(',idari(n),')=',snm(2),'  ,M(',idari(n),')=',snm(3)
    print'(A,A,I2,A,f15.10,A,I2,A,f15.10,A,I2,A,f15.10)','             ',&
        &'   S(',migi(n),')=',snm(4),'  ,N(',migi(n),')=',snm(5),'  ,M(',migi(n),')=',snm(6)
   !
   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(6,6), b(6,6), c(6,6)
!
!
do i=1,6
do j=1,6
c(i,j)=0.
 do k=1,6
  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(6,6)
do i=1,3
do j=1,3
 t(i+3,j)=0.d0
 t(i,j+3)=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)=1.d0
t(4,4)=t(1,1)
t(4,5)=t(1,2)
t(5,4)=t(2,1)
t(5,5)=t(2,2)
t(6,6)=1.d0
return
end
!
subroutine mxtmx(a,b,c)
implicit real*8(a-h, o-z)
dimension a(6,6), b(6,6), c(6,6)
!
do i=1,6
do j=1,6
c(i,j)=0.
 do k=1,6
  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(45,45),x(45),b(45)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす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
!
!

