implicit real*8(a-h, o-z)
dimension f1(4),s1(4,4),d1(4),f2(4),s2(4,4),d2(4),t(4,4),&
& f(6),s(6,6),d(6),a1(4,4),a2(4,4),x(6)
do i=1,4
do j=1,4
s1(i,j)=0.d0
s2(i,j)=0.d0
 t(i,j)=0.d0
end do
end do
! 要素1の剛性マトリクス
s1(2,2)= 1.d0; s1(2,4)=-1.d0
s1(4,2)=-1.d0; s1(4,4)= 1.d0
sk2=1.d0/sqrt(2.d0)
! 要素2の剛性マトリクス
s2(2,2)= sk2; s2(2,4)=-sk2
s2(4,2)=-sk2; s2(4,4)=sk2
! 座標変換マトリクス
t(1,1)=-sk2; t(1,2)=-sk2
t(2,1)= sk2; t(2,2)=-sk2
t(3,3)=-sk2; t(3,4)=-sk2
t(4,3)= sk2; t(4,4)=-sk2
!
do i=1,6
x(i)=1.d0 !境界条件に1を入れてしまう
f(i)=0.d0
end do
! 境界条件をxに入れる
! 拘束する節点に0.d0を入れる
x(1)=0.d0; x(2)=0.d0
x(5)=0.d0; x(6)=0.d0
! 荷重条件
f(3)=1.d0
!
!print*,'s1='
!do i=1,4
!print'(4f10.3)', (s1(i,j),j=1,4)
!end do
!!
!print*,'s2='
!do i=1,4
!print'(4f10.3)', (s2(i,j),j=1,4)
!end do
!
call mxtmx(t,s2,a1)
!!
!print*,'a1='
!do i=1,4
!print'(4f10.3)', (a1(i,j),j=1,4)
!end do
!
call mxmx(a1,t,a2)
!!
!print*,'tk2t='
!do i=1,4
!print'(4f10.3)', (a2(i,j),j=1,4)
!end do
!
!
do i=1,6
do j=1,6
s(i,j)=0.d0
end do
end do
!
do i=1,4
do j=1,4
s(i,j)=s1(i,j)
end do
end do
!
do i=1,4
do j=1,4
s(i+2,j+2)=s(i+2,j+2)+a2(i,j)
end do
end do
!
!!
!　境界条件を入れる
do i=1,6
do j=1,6
  s(i,j)=x(i)*s(i,j)
  s(j,i)=x(i)*s(j,i)
end do
end do
do i=1,6
if(x(i)<1.d-3) then
    s(i,i)=1.d0
end if
end do
!!
!!
!!
print*,'s='
do i=1,6
print'(6f10.3)', (s(i,j),j=1,6)
end do
!
!
call gausu(6,s,d,f)
print*,'d='
print'(6f10.3)', (d(j),j=1,6)
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)
!
!!
print*,'a1='
do i=1,4
print'(4f10.3)', (a(i,j),j=1,4)
end do
!
!!
print*,'t='
do i=1,4
print'(4f10.3)', (b(i,j),j=1,4)
end do
!
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
!
!!
print*,'c='
do i=1,4
print'(4f10.3)', (c(i,j),j=1,4)
end do
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)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15
!
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
!
!

