主程序main.f90,搞了两个例子进行调用,例子是钟院士那篇文献里的,附后的文献里有的。一个程序是reciseTIM.f90,里面包含了精细时程积分的子程序。
最后的PDF文件是程序说明.
主程序Program Main
- use PreciseTimeIntegration
- Real(8), Dimension(8, 8) :: K,M,C,INVK
- Real(8), Dimension(8, 4) :: R,U,V,A
- Real(8) dt
- dt=1
- R=0
- C=0
- C(7,7)=5
- C(7,8)=-5
- C(8,7)=-5
- C(8,8)=5
- DO i=1,8
- M(i,i)=8
- END DO
- DO i=1,7
- K(i,i+1)=-4
- k(i+1,i)=-4
- END DO
- DO i=2,7
- k(i,i)=8
- END DO
- K(1,1)=4
- K(8,8)=4
- ! WRITE (*,*) M
- ! WRITE (*,*) K
- ! WRITE (*,*) C
- ! INVK=INV(K,8)
- ! STOP
- U=0
- U(8,1)=10
- V=0
- A=0
- call PreciseTIM (M, C, k, R, 8, 4, dt, U, V, A)
- write (*,*) U
- ! write (*,*) V
- ! write (*,*) A
- end
复制代码
本程序来自http://bbs.tongji.edu.cn/bbsanc.php?path=%2FPersonalCorpus%2FA%2FABAYA%2FD4F47E4DC%2FA41E7F755
精细时程积分文件:
PreciseTIM.f90
!-----------------------------------------------------------------
Module PreciseTimeIntegration
Contains
- !//////////////////////////////////////////////////////////////
- Subroutine PreciseTIM (M, G, K, R, N, Nt, tao, X, XX, XXX)
- !//////////////////////////////////////////////////////////////
- ! 精细时程积分法程序,
- ! ABAYA,同济大学建筑工程系,2001.7.20
- !!//////////////////////////////////////////////////////////////
- Real*8 M(N,N), G(N,N), K(N,N), R(N,Nt), X(N,Nt), &
- XX(N,Nt), XXX(N,Nt), H(2*N,2*N), T(2*N,2*N), &
- R0(2*N,Nt), R1(2*N,Nt), B(N,N), C(N,N), tao
- Call CalH (M, G, K, N, Nt, H, B, C)
- Call CalT (H, tao, N, T)
- Call CalR0R1 (R, R0, R1, N, Nt)
- Call CalX (T, H, R, R0, R1, tao, N, Nt, M, G, B, C, X, XX, XXX)
- End Subroutine PreciseTIM
- !//////////////////////////////////////////////////////////////
- Subroutine CalH (M, G, K, N, Nt, H, B, C)
- Real*8 M(N,N), G(N,N), K(N,N), R(N,Nt), INVM(N,N), &
- A(N,N), B(N,N), C(N,N), D(N,N), H(2*N,2*N)
- INVM=INV(M,N)
- A=-0.5*Matmul(INVM,G)
- B=0.25*Matmul(Matmul(G,INVM),G)-K
- C=-0.5*Matmul(G,INVM)
- D=INVM
- H(1:N,1:N)=A
- H(N+1:2*N,1:N)=B
- H(1:N,N+1:2*N)=D
- H(N+1:2*N,N+1:2*N)=C
- End Subroutine CalH
- !//////////////////////////////////////////////////////////////
- Subroutine CalT (H, tao, N, T)
- Real*8, Dimension(2*N,2*N) :: H, T, Ta, I
- Real*8 tao, dt
- Integer m
- m=2**20
- dt=tao/m
- I=0; T=0
- Do j=1,2*N
- I(j,j)=1
- End do
- Ta=Matmul(dt*H,(I+0.5*dt*H))
- Do j=1,20
- Ta=2*Ta+Matmul(Ta,Ta)
- End do
- T=I+Ta
- End Subroutine CalT
- !//////////////////////////////////////////////////////////////
- Subroutine CalR0R1 (R, R0, R1, N, Nt)
- Real*8 R(N,Nt), R0(2*N,Nt), R1(2*N,Nt)
- R0=0; R1=0;
- R0(N+1:2*N,:)=R
- Do i=2,Nt
- R1(N+1:2*N,i-1)=R(:,i)-R(:,i-1)
- End do
- End Subroutine CalR0R1
- !//////////////////////////////////////////////////////////////
- Subroutine CalX (T, H, R, R0, R1, tao, N, Nt, M, G, B, C, X, XX, XXX)
- Real*8 X(N,Nt), XX(N,Nt), XXX(N,Nt), H(2*N,2*N), &
- T(2*N,2*N), R0(2*N,Nt), R1(2*N,Nt), V(2*N,Nt), &
- p(N,Nt), q(N,Nt), R(N,Nt), B(N,N), C(N,N), M(N,N), &
- G(N,N), tao, INVH(2*N,2*N), index
-
- q=X
- p=Matmul(M,XX)+0.5*Matmul(G,X)
- V(1:N,:)=q
- V(N+1:2*N,:)=p
- index=0
- Do i=1,2*N
- Do j=1,Nt
- IF (ABS(R0(i,j)).GT.1E-8) THEN
- index=1; GOTO 10
- END IF
- End do
- End do
- 10 If (abs(index).gt.1e-8) then
- INVH=INV(H,2*N)
- Do i=2,Nt
- WRITE (*,'("*********** LOAD STEP: ",I5," ***********")') i
- V(:,i)=Matmul(T,(V(:,i-1)+Matmul(INVH,(R0(:,i-1)+ &
- Matmul(INVH,R1(:,i-1))))))-Matmul(INVH &
- ,(R0(:,i-1)+Matmul(INVH,R1(:,i-1))+ &
- tao*R1(:,i-1)))
- End do
- Else
- Do i=2,Nt
- WRITE (*,'("*********** LOAD STEP: ",I5," ***********")') i
- V(:,i)=Matmul(T,V(:,i-1))
- End do
- End if
- q=V(1:N,:)
- p=V(N+1:2*N,:)
- Do i=1,Nt
- X(:,i)=q(:,i)
- XX(:,i)=Matmul(INV(M,N),p(:,i))-0.5*Matmul(Matmul(INV(M,N),G),X(:,i))
- XXX(:,i)=Matmul(Matmul(INV(M,N),B),X(:,i))-0.5*Matmul(Matmul(INV(M,N) &
- ,G),XX(:,i))+Matmul(Matmul(INV(M,N),C),p(:,i))+Matmul(INV(M,N),R(:,i))
- End do
- End Subroutine CalX
- !//////////////////////////////////////////////////////////////
- Function INV (A, N)
- Use Numerical_Libraries
- REAL*8, DIMENSION(N,N) :: A, INV
- INV=0
- CALL DLINRG (N, A, N, INV, N)
- End Function INV
- !//////////////////////////////////////////////////////////////
复制代码
End Module PreciseTimeIntegration |