! 文字コードはutf-8
! 自由形式のfortranで書いています。
! g77でコンパイルする場合は、
!g77 -ffree-form -o smv smv.f
! のようにコンパイルしてください。
! gfortran でコンパイルする場合は文字コードをutf-8にしてください。
program smv
implicit real*8(a-h,o-z)
dimension s122122(122,122),d122(122),&  !全体剛性行列
&    f4(4),s44(4,4),sk44(4,4),d(4),& !要素剛性行列
& bs61(61),bm61(61),b122(122),&  !境界条件ベクトル
& fs61(61),fm61(61),f122(122) !荷重条件ベクトル
!
!
ell=1.d0 ! 梁の長さL
ei=1.d0  !曲げ剛性EI
!
!
do i=1,61
bs61(i)=1.d0
bm61(i)=1.d0
fs61(i)=0.d0
fm61(i)=0.d0
end do
! 
! 以下から境界条件を選択または自分で設定 
!call tanzyunb(bs61,bm61) !単純支持     △------------△
!call nitten3b(bs61,bm61) !1/3点2箇所支持 ----△------△---
!call satankotei(bs61,bm61) !左端固定    |------------
!call utankotei(bs61,bm61) !右端固定    ------------|
!call ryoukotei(bs61,bm61) !両端固定    |------------|
call sakouhin(bs61,bm61) !左端固定右端ヒンジ  |------------△
!call sahinuko(bs61,bm61) !左端ヒンジ右端固定  △------------|
!call sansizi(bs61,bm61) ! 3点支持     △------△------△
!call sansizi43(bs61,bm61) ! 3点支持     △---------△---△
!call haridasi(bs61,bm61) ! 2点支持     △------△------
!call haridasig(bs61,bm61) ! 2点支持(逆対称4点曲げ用) △---3L---△-L-
! 
! 以下から荷重条件を選択または自分で設定 
!call toubunpuf(fs61,fm61) !等分布荷重 ↓↓↓↓↓↓↓↓↓
!call toubunpuuf(fs61,fm61) !等分布荷重上向き ↑↑↑↑↑↑↑↑↑
!call hansaif(fs61,fm61) !半載荷重 ↓↓↓↓--------
!call hansaif43(fs61,fm61) !半載荷重 ↓↓↓↓--
!call hansaifm(fs61,fm61) !右側半載荷重 ------↓↓↓↓↓
!call sankakuf(fs61,fm61) !三角形分布荷重 0____----^^^^1↓
!call sankakusf(fs61,fm61) !三角形分布荷重右下がり 1^^^^----____0↓
!call sankakucf(fs61,fm61) !中央三角形分布荷重 ____0_-^1↓____ 
call tyuuouf(fs61,fm61) !中央鉛直荷重載荷  ------↓------
!call tyuuoufuu(fs61,fm61) !中央鉛直右端上向き載荷   ------↓-----↑
!call tyuuouuf(fs61,fm61) !中央上向き鉛直荷重載荷  ------↑------
!call tyuuoumm(fs61,fm61) !中央右回りモーメント  ------×------
!call tyuuoumh(fs61,fm61) !中央左回りモーメント  ------◎------
!call hidarimm(fs61,fm61) !左端右回りモーメント  ×------------
!call migimm(fs61,fm61) !右端右回りモーメント  ------------×
!call ryoutanf(fs61,fm61) !両端載荷 ↓-----------------↓
!call ryoutanuf(fs61,fm61) !両端載荷上向き ↑-----------------↑
!call ryoutanmf(fs61,fm61) !両端載荷右回り ↑----------------↓
!call itten3f(fs61,fm61) !1/3点1箇所載荷  ---↓---------
!call nitten3f(fs61,fm61) !1/3点2箇所載荷  ----↓----↓----
!call nitten3fh(fs61,fm61) !1/3点2箇所載荷左回り  ----↓----↑----
!call niten4f(fs61,fm61) ! 1/4点2箇所載荷   ---↓------↓---
!call niten4hf(fs61,fm61) ! 1/4点2箇所載荷張出用   ---↓2--------↓1
!call santen4f(fs61,fm61) ! 1/4点3箇所載荷   ---↓---↓---↓---
!call yonten5f(fs61,fm61) ! 1/5点4箇所載荷   --↓--↓--↓--↓--
! 
call sm61sm122(bs61,bm61,b122,fs61,fm61,f122) 
!
!
! 1要素の要素剛性行列 s44
s44(1,1)=12.d0;s44(1,2)=-6.d0;s44(1,3)=-s44(1,1);s44(1,4)=s44(1,2)
               s44(2,2)= 4.d0;s44(2,3)=-s44(1,2);s44(2,4)=2.d0
                              s44(3,3)= s44(1,1);s44(3,4)=s44(2,3)
                                                 s44(4,4)=s44(2,2)
!重ね合わせられる要素剛性行列 sk44
 do i=1,2
 do j=i,4
   sk44(i,j)=s44(i,j)
 end do
 end do
!
 sk44(1,1)=s44(1,1)+s44(3,3)
 sk44(1,2)=s44(1,2)+s44(3,4)
 sk44(2,2)=s44(2,2)+s44(4,4)
!
 do i=2,4
 do j=1,i-1
 s44(i,j)=s44(j,i)
 sk44(i,j)=sk44(j,i)
 end do
 end do
!
!
call zentai(s44,sk44,s122122,b122)
!do i=1,4
!print*,(sk44(i,j),j=1,4)
!end do
!do i=1,10
!print*,(s122122(i,j),j=1,10)
!end do
!print*
!do i=113,122
!print*,(s122122(i,j),j=113,122)
!end do
!
call gausu(122,s122122,d122,f122)
!print*,((s44(i,j),j=1,4),'\n',i=1,4)
!
! たわみの出力
open(7,file='tawami.d')
do i=1,121,2
write(7,*) (real(i)+1.d0)/2.d0,d122(i)
end do
close(7)
!
! 軸線の出力 
open(7,file='ziku.d')
write(7,*)  1,0
write(7,*) 61,0
close(7)
!
call nairyoku(s44,d122)
!
end
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
subroutine zentai(s44,sk44,s122122,b122)
implicit real*8(a-h,o-z)
dimension s44(4,4),sk44(4,4),s122122(122,122),b122(122)
!
 do i=1,122
 do j=1,122
    s122122(i,j)=0.d0
 end do
 end do
! 1要素目の2行目まで
do i=1,2
do j=i,4
   s122122(i,j)=s44(i,j)
end do
end do
!
! 1要素目の3行目から60要素目の2行目まで
do ie=1,59
do i=1,2
do j=i,4
   s122122(2*ie+i,2*ie+j)=sk44(i,j)
end do
end do
end do
!
! 60要素目の3行目から
do i=3,4
do j=i,4
   s122122(118+i,118+j)=s44(i,j)
end do
end do
!
!
do i=2,122
do j=1,i-1
 s122122(i,j)=s122122(j,i)
end do
end do
!
!
do i=1,122
do j=1,122
   s122122(i,j)=b122(i)*s122122(i,j)
   s122122(j,i)=b122(i)*s122122(j,i)
end do
end do
!
do i=1,122
   if( b122(i)<1.d-3 ) then
           s122122(i,i)=1.d0
   end if
end do
!
!
return
end
!
!
subroutine sm61sm122(bs61,bm61,b122,fs61,fm61,f122) 
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61),b122(122),fs61(61),fm61(61),f122(122)
do i=1,61
b122(2*i-1)=bs61(i)
b122(2*i  )=bm61(i)
!
f122(2*i-1)=fs61(i)
f122(2*i  )=fm61(i)
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
!
!
! 内力(せん断力、曲げモーメント)の計算
subroutine nairyoku(s44,d122)
implicit real*8(a-h,o-z)
dimension s44(4,4),d122(122),fn122(122)
!
do i=1,122
fn122(i)=0.d0
end do
!
! 1から60要素までの左端の内力
do n=1,60
do i=1,2
do j=1,4
   fn122(n*2-2+i)=fn122(n*2-2+i)-s44(i,j)*d122(n*2-2+j)
end do
end do
end do
!
! 60要素目の右端の内力
do i=3,4
do j=1,4
   fn122(118+i)=fn122(118+i)+s44(i,j)*d122(118+j)
end do
end do
!
open(7,file='sendan.d')
write(7,*) 1,0.
do i=1,60
write(7,*) i,fn122(i*2-1)
end do
write(7,*) 61,fn122(121)
write(7,*) 61,0.
close(7)
!
!
open(7,file='sendant.d')
do i=1,60
write(7,*) i,fn122(i*2-1)
end do
write(7,*) 61,fn122(121)
close(7)
!
open(7,file='mage.d')
do i=1,60
write(7,*) i,fn122(i*2)
end do
write(7,*) 61,fn122(122)
close(7)
!
!
return
end
!
!
!
!
! 単純支持
subroutine tanzyunb(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(61)=0.d0
return
end
!
! 左端固定
subroutine satankotei(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bm61(1)=0.d0
return
end
!
! 右端固定
subroutine utankotei(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(61)=0.d0
bm61(61)=0.d0
return
end
!
! 両端固定
subroutine ryoukotei(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bm61(1)=0.d0
bs61(61)=0.d0
bm61(61)=0.d0
return
end
!
! 左端固定右端ヒンジ
subroutine sakouhin(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bm61(1)=0.d0
bs61(61)=0.d0
return
end
!
! 左端ヒンジ右端固定
subroutine sahinuko(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(61)=0.d0
bm61(61)=0.d0
return
end
!
! 3点固定
subroutine sansizi(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(31)=0.d0
bs61(61)=0.d0
return
end
!
!
! 3点固定
subroutine sansizi43(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(43)=0.d0
bs61(61)=0.d0
return
end
!
!
! 2点固定張り出し
subroutine haridasi(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(31)=0.d0
return
end
!
! 2点固定張り出し(逆対称4点曲げ用)
subroutine haridasig(bs61,bm61)
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(1)=0.d0
bs61(46)=0.d0
return
end
!
subroutine nitten3b(bs61,bm61) !1/3点2箇所支持 ----△------△---
implicit real*8(a-h,o-z)
dimension bs61(61),bm61(61)
bs61(21)=0.d0
bs61(41)=0.d0
return
end
!
!   鉛直荷重外力はfs61にPL**2/EIの次元で与える
! モーメント外力はfm61に   PL/EIの次元で与える
!
! 等分布荷重
subroutine toubunpuf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,61
fs61(i)=1.d0
end do
return
end
!
!
! 等分布荷重上向き
subroutine toubunpuuf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,61
fs61(i)=-1.d0
end do
return
end
!
!
! 半載荷重
subroutine hansaif(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,31
fs61(i)=1.d0
end do
return
end
! 半載荷重
subroutine hansaif43(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,43
fs61(i)=1.d0
end do
return
end
!
!
!右側半載荷重 ------↓↓↓↓↓
subroutine hansaifm(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=31,61
fs61(i)=1.d0
end do
return
end
!
! 三角形分布荷重 0____----^^^^1↓
subroutine sankakuf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,61
fs61(i)=real(i)/61
end do
return
end
!
! 三角形分布荷重右下がり 1^^^^----____0↓
subroutine sankakusf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=1,61
fs61(i)=1.d0-real(i-1)/61
end do
return
end
!
subroutine sankakucf(fs61,fm61) !中央三角形分布荷重 ____0_-^1↓____ 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
do i=21,41
fs61(i)=(real(i)-21.d0)/20
end do
return
end
!
! 中央1点載荷
subroutine tyuuouf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(31)=1.d0
return
end
!
subroutine tyuuoufuu(fs61,fm61) !中央鉛直右端上向き載荷   ------↓-----↑
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(31)=1.d0
fs61(61)=-1.d0
return
end
!
! 中央1点載荷(上向き)
subroutine tyuuouuf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(31)=-1.d0
return
end
!
!中央右回りモーメント  ------×------
subroutine tyuuoumm(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fm61(31)=-1.d0
return
end
!
!中央左回りモーメント  ------◎------
subroutine tyuuoumh(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fm61(31)=1.d0
return
end
!
!左端右回りモーメント  ×------------
subroutine hidarimm(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fm61(1)=-1.d0
return
end
!
!
!右端右回りモーメント  ------------×
subroutine migimm(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fm61(61)=-1.d0
return
end
!
!
subroutine ryoutanf(fs61,fm61) !両端載荷 ↓-----------------↓
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(1)=1.d0
fs61(61)=1.d0
return
end
!
subroutine ryoutanuf(fs61,fm61) !両端載荷上向き ↑-----------------↑
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(1)=-1.d0
fs61(61)=-1.d0
return
end
!
subroutine ryoutanmf(fs61,fm61) !両端載荷右回り ↑----------------↓
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(1)=-1.d0
fs61(61)=1.d0
return
end
!
!1/3点1箇所載荷  ---↓---------
subroutine itten3f(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(21)=1.d0
return
end
!
!1/3点2箇所載荷  ----↓----↓----
subroutine nitten3f(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(21)=1.d0
fs61(41)=1.d0
return
end
!
!1/3点2箇所載荷左回り  ----↓----↑----
subroutine nitten3fh(fs61,fm61) 
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(21)=1.d0
fs61(41)=-1.d0
return
end
!
! 1/4点2箇所載荷 
subroutine niten4f(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(16)=1.d0
fs61(46)=1.d0
return
end
!
!
! 1/4点2箇所載荷 張り出し梁用
subroutine niten4hf(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(16)=2.d0
fs61(61)=1.d0
return
end
!
! 1/4点3箇所載荷 
subroutine santen4f(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(16)=1.d0
fs61(31)=1.d0
fs61(46)=1.d0
return
end
!
! 1/5点4箇所載荷 
subroutine yonten5f(fs61,fm61)
implicit real*8(a-h,o-z)
dimension fs61(61),fm61(61)
fs61(13)=1.d0
fs61(25)=1.d0
fs61(37)=1.d0
fs61(49)=1.d0
return
end
!
