[分享]精细时程积分FORTRAN源程序
主程序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
10If (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 赞一个先,
自己也曾编过线性的wilson时程分析程序.
可以解决剪切层模型问题.没有作验证:) 对此比较感兴趣,先睹为快 谢谢了,感谢共享这么好的源代码。 谢谢共享 pdf的文档怎么无法下载? 使用marc进行有限元分析时一定要用fortran语言进行有关的编程吗?可不可以用C语言 Marc的接口是Fortran的,会C语言学Fortran不难
Marc各个版本与fortran版本搭配如下:
marc2003需要fortran6.0以上版本
marc2000以及marc2001需要fortran5.0以上版本 ding 哈哈,找拉好久的拉
好!!!! 非常不错的东西,感谢搂住的热心 好咚咚 先看一下,是否有用 xiexie