📄 eq.for
字号:
END DO
VV(G_ELE(I,1)*2) = VV(G_ELE(I,1)*2)-T1
VV(G_ELE(I,2)*2) = VV(G_ELE(I,2)*2)-T2
VV(G_ELE(I,3)*2) = VV(G_ELE(I,3)*2)-T3
VV(G_ELE(I,4)*2) = VV(G_ELE(I,4)*2)-T4
DO P = 1,MTYPE
IF (DEVIDE(P) .EQ. 1) THEN
VV(G_ELE(ELE,4+P)*2)=VV(G_ELE(ELE,4+P)*2)-TEMP(P)
END IF
END DO
ELSE IF (MTYPE .EQ. 6) THEN
IF (N_INT .EQ. 1) THEN
TT=1
ELSE IF(N_INT .EQ. 2) THEN
TT = 3
ELSE IF(N_INT .EQ. 3) THEN
TT = 4
ELSE
PRINT*,'ERROR N_INT IN SUBROUTINE G_LOAD'
END IF
CALL CAL_POINT()
DO K = 1,TT
KSI = POINT3(K,1)
YITA = POINT3(K,2)
CALL CAL_N()
DO P = 1,MTYPE
VV(G_ELE(I,P)*2)=VV(G_ELE(I,P)*2)-
5 H3(K)*ROU*GG*PT*N(P)
END DO
END DO
ELSE IF (MTYPE .EQ. 8) THEN
CALL CAL_POINT()
DO J = 1,N_INT
DO K = 1,N_INT
YITA = POINT(J)
KSI = POINT(K)
CALL CAL_N()
DO P = 1,MTYPE
VV(G_ELE(I,P)*2) = VV(G_ELE(I,P)*2)
4 -H(K)*H(J)*PT*ROU*GG*N(P)
END DO
END DO
END DO
END IF
END DO
END SUBROUTINE G_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE EQUAL_VERTICAL_LINE_LOAD()
!CALCULATE THE LOAD (VERTICAL LOAD ON LINE EQUALLY)
USE DEF
IMPLICIT NONE
!ABOVE TO TRANSFER
INTEGER(4)::I,J,T(1:2),GI,GJ
REAL(8)::Q,YI,YJ,XI,XJ
DO I = 1,N_LINE_LOAD
Q = LINE_LOAD(I)
T(1) = LINE_LOAD_NODE(I,1)
T(2) = LINE_LOAD_NODE(I,2)
XI = G_NODE(T(1),1)
XJ = G_NODE(T(2),1)
YI = G_NODE(T(1),2)
YJ = G_NODE(T(2),2)
DO J = 1,2
GI = 2*T(J)-1
GJ = 2*T(J)
VV(GI) = VV(GI)+Q*PT*(YI-YJ)/2.0
VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)/2.0
END DO
END DO
END SUBROUTINE EQUAL_VERTICAL_LINE_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE EQUAL_X_LOAD()
!CALCULATE THE LOAD OF KIND
USE DEF
IMPLICIT NONE
!ABOVE TO TRANSFER
INTEGER(4)::T1,T2,I
REAL(8)::XI,XJ,YI,YJ,DIST,Q
DO I = 1,X_N_LOAD
Q = X_LOAD(I)
T1 = X_LOAD_NODE(I,1)
T2 = X_LOAD_NODE(I,2)
XI = G_NODE(T1,1)
XJ = G_NODE(T2,1)
YI = G_NODE(T1,2)
YJ = G_NODE(T2,2)
DIST = SQRT((XI-XJ)**2+(YI-YJ)**2)
VV(2*T1-1) = VV(2*T1-1)+Q*PT*DIST/2.0
VV(2*T2-1) = VV(2*T2-1)+Q*PT*DIST/2.0
END DO
END SUBROUTINE EQUAL_X_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE LINEAR_LOAD()
!CALCULATE THE LINEARLY LOAD ON ONE VERTICLE LINE
USE DEF
IMPLICIT NONE
!ABOVE TO TRANSFER
INTEGER(4)::I,J,NN,T(1:2),GI,GJ
REAL(8)::Q,XI,XJ,YI,YJ
!N_LINEAR(1:3):1.MINIMAL NUMBER OF STARTING POINT WITH THE LINEAR FORCE
! 2.MAXIMUM NUMBER OF ENDDING POINT WITH THE LINEAR FORCE
!Q_LINEAR:THE MAXIMUM VALUE OF THE FORCE
! Q_LINEAR .GT. 0: PRESSURE
! .LT. 0: PULL
Q = Q_LINEAR
NN = N_LINEAR(2)-N_LINEAR(1)
DO I = 1,NN
T(1) = N_LINEAR(1)+(I-1)
T(2) = T(1)+1
XI = G_NODE(T(1),1)
XJ = G_NODE(T(2),1)
YI = G_NODE(T(1),2)
YJ = G_NODE(T(2),2)
GI = 2*T(1)-1
GJ = 2*T(1)
VV(GI) = VV(GI)+Q*PT*(YI-YJ)*(1/3.0+(I-1)/NN/2.0)
VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)*(1/3.0+(I-1)/NN/2.0)
GI = 2*T(2)-1
GJ = 2*T(2)
VV(GI) = VV(GI)+Q*PT*(YI-YJ)*(1/6.0+(I-1)/NN/2.0)
VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)*(1/6.0+(I-1)/NN/2.0)
END DO
END SUBROUTINE LINEAR_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ADD_BC()
!ADD THE BOUNDARY CONDITION
USE DEF
IMPLICIT NONE
!ABOVE TO TRANSFER
INTEGER(4)::I,J,GI,GJ
REAL(8)::ALPHA
ALPHA = 1.0E10
DO I = 1,N_BC
GI = 2*N_DIS(I)-1
GJ = 2*N_DIS(I)
IF (BC(I,1) .NE. 1.0) THEN
KK(DIAG(GI)) = ALPHA*KK(DIAG(GI))
VV(GI) = ALPHA*BC(I,1)*KK(DIAG(GI))
END IF
IF (BC(I,2) .NE. 1.0) THEN
KK(DIAG(GJ)) = ALPHA*KK(DIAG(GJ))
VV(GJ) = ALPHA*BC(I,2)*KK(DIAG(GJ))
END IF
END DO
END SUBROUTINE ADD_BC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SOLVE(A,P,M,D,NUM)
!一维变带宽矩阵求解方程组
USE DEF
INTEGER(4)::NUM
REAL(8)::A(DOF)
REAL(8)::D(NUM)
INTEGER(4)::P(NUM),M(NUM+1)
!以上变量需要传入
INTEGER(4)::I,J,Q
REAL(8)::S,T
!用一维数组A存储上三角的非零元素
!用一维数组m储存主对角元在一维数组中的位置
!用一维数组N储存列高
!用一维数组p储存每一列的起始行号
!GUASS分解a
IF (1 .EQ. 0) THEN
DO J=1,NUM
DO I = 1,J
S = A(M(J)+J-I)
DO Q = MAX(P(I),P(J)), I-1
S=S-A(M(I)+I-Q)*A(M(J)+J-Q)/A(M(Q))
END DO
A(M(J)+J-I) = S
END DO
END DO
!求解方程
!D是最后求出来的解
DO J =1,NUM
S = D(J)
DO I = P(J),J-1
S = S-A(M(J)+J-I)*D(I)/A(M(I))
END DO
D(J) = S
END DO
D(NUM) = D(NUM)/A(M(NUM))
DO I = NUM-1,1,-1
S = D(I)
DO J = I+1,NUM
IF (P(J) .LE. I) THEN
S = S-A(M(J)+J-I)*D(J)
END IF
END DO
D(I) = S/A(M(I))
END DO
END IF
!LDLT
IF (1 .EQ. 1) THEN
DO K = 1,NUM
DO J = P(K),K-1
A(M(K)) = A(M(K)) - A(M(K)+K-J)**2*A(M(J))
END DO
DO I = K+1,NUM
IF (P(I) .LE .K) THEN
DO J = MAX(P(K),P(I)),K-1
A(M(I)+I-K)=A(M(I)+I-K)-A(M(I)+I-J)*A(M(K)+K-J)*A(M(J))
END DO
A(M(I)+I-K)=A(M(I)+I-K)/A(M(K))
END IF
END DO
END DO
DO K = 2,NUM
DO J = P(K),K-1
D(K) = D(K) - A(M(K)+K-J)*D(J)
END DO
END DO
D(NUM)=D(NUM)/A(M(NUM))
DO K = NUM-1,1,-1
T = 0.0D0
DO J = K+1,NUM
IF (P(J) .LE. K) THEN
T = T+A(M(J)+J-K)*D(J)
END IF
END DO
D(K) = D(K)/A(M(K))-T
END DO
END IF
DO I = 1,NUM
IF (ABS(D(I)) .LT. 1.0E-16) THEN
D(I) = 0.0
END IF
END DO
END SUBROUTINE SOLVE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE CAL_STR()
!CALCULATE THE STRESS
USE DEF
IMPLICIT NONE
!ABOVE TO TRANSFER
INTEGER(4)::I,J,K,P,Q,MID
REAL(8)::DD(3,3),BB(3,16)
REAL(8)::QQ1,QQ2,QQ3,TEMP1(1:16),TEMP2(1:3)
EXTERNAL::MULTI,CAL_BB,CAL_DD
SIGMA=0.0
DO I = 1,N_ELE
BB = 0.0
DD = 0.0
MTYPE = MARK(I)
ELE = I
COORD = 0.0
MID = N_INT
N_INT = 2
QQ1 = 0.0D0
QQ2 = 0.0D0
QQ3 = 0.0D0
DO J = 1,8
IF (G_ELE(I,J) .NE. 0) THEN
COORD(J,1) = G_NODE(G_ELE(I,J),1)
COORD(J,2) = G_NODE(G_ELE(I,J),2)
END IF
END DO
IF (MTYPE .EQ. 3) THEN
CALL FIND_34()
CALL CAL_POINT()
DO K = 1,MTYPE
KSI = POINT3(K,1)
YITA = POINT3(K,2)
CALL CAL_DN()
CALL COMPUT()
DO J = 1,MTYPE
IF (DEVIDE(J) .EQ. 1) THEN
DN(1,J) = DN6(1,J)
DN(2,J) = DN6(2,J)
DN(1,3+J) = DN6(1,3+J)
DN(2,3+J) = DN6(2,3+J)
END IF
END DO
CALL CAL_BB(BB)
CALL CAL_DD(DD)
CALL MULTI(DD,BB,BB,3,3,16)
TEMP1 = 0.0
DO J = 1,16
IF (G_ELE(I,J) .NE. 0) THEN
TEMP1(2*J-1) = VV(2*G_ELE(I,J)-1)
TEMP1(2*J) = VV(2*G_ELE(I,J))
END IF
END DO
CALL MULTI(BB,TEMP1,TEMP2,3,16,1)
QQ1 = QQ1+TEMP2(1)
QQ2 = QQ2+TEMP2(2)
QQ3 = QQ3+TEMP2(3)
END DO
ELSE IF (MTYPE .EQ. 4) THEN
CALL FIND_34()
CALL CAL_POINT()
DO K = 1,2
DO P = 1,2
KSI = POINT(K)
YITA = POINT(P)
CALL CAL_DN()
CALL COMPUT()
DO J = 1,MTYPE
IF (DEVIDE(J) .EQ. 1) THEN
DN(1,J) = DN6(1,J)
DN(2,J) = DN6(2,J)
DN(1,3+J) = DN6(1,3+J)
DN(2,3+J) = DN6(2,3+J)
END IF
END DO
CALL CAL_BB(BB)
CALL CAL_DD(DD)
CALL MULTI(DD,BB,BB,3,3,16)
TEMP1 = 0.0
DO J = 1,16
IF (G_ELE(I,J) .NE. 0) THEN
TEMP1(2*J-1) = VV(2*G_ELE(I,J)-1)
TEMP1(2*J) = VV(2*G_ELE(I,J))
END IF
END DO
CALL MULTI(BB,TEMP1,TEMP2,3,16,1)
QQ1 = QQ1+TEMP2(1)
QQ2 = QQ2+TEMP2(2)
QQ3 = QQ3+TEMP2(3)
END DO
END DO
ELSE IF (MTYPE .EQ. 6) THEN
CALL CAL_POINT()
DO K = 1,3
KSI = POINT3(K,1)
YITA = POINT3(K,2)
CALL CAL_DN()
CALL COMPUT()
CALL CAL_BB(BB)
CALL CAL_DD(DD)
CALL MULTI(DD,BB,BB,3,3,16)
TEMP1 = 0.0
DO J = 1,8
IF (G_ELE(I,J) .NE. 0) THEN
TEMP1(2*J-1) = VV(2*G_ELE(I,J)-1)
TEMP1(2*J) = VV(2*G_ELE(I,J))
END IF
END DO
CALL MULTI(BB,TEMP1,TEMP2,3,16,1)
QQ1 = QQ1+TEMP2(1)
QQ2 = QQ2+TEMP2(2)
QQ3 = QQ3+TEMP2(3)
END DO
ELSE IF (MTYPE .EQ. 8) THEN
CALL CAL_POINT()
DO K = 1,2
DO P = 1,2
KSI = POINT(K)
YITA = POINT(P)
CALL CAL_DN()
CALL COMPUT()
CALL CAL_BB(BB)
CALL CAL_DD(DD)
CALL MULTI(DD,BB,BB,3,3,16)
TEMP1 = 0.0
DO J = 1,8
IF (G_ELE(I,J) .NE. 0) THEN
TEMP1(2*J-1) = VV(2*G_ELE(I,J)-1)
TEMP1(2*J) = VV(2*G_ELE(I,J))
END IF
END DO
CALL MULTI(BB,TEMP1,TEMP2,3,16,1)
QQ1 = QQ1+TEMP2(1)
QQ2 = QQ2+TEMP2(2)
QQ3 = QQ3+TEMP2(3)
END DO
END DO
ELSE
PRINT*,'WRONG ELEMENT TYPE IN SUBROUTINE CAL_STR'
END IF
SIGMA(ELE,1)=QQ1/3.0D0
SIGMA(ELE,2)=QQ2/3.0D0
SIGMA(ELE,3)=QQ3/3.0D0
N_INT = MID
END DO
END SUBROUTINE CAL_STR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM MAIN
USE DEF
IMPLICIT NONE
INTEGER(4)::I,J,K
EXTERNAL::READ_DATA,CAL_PE,BAND_K,FORM_LOAD,ADD_BC,SOLVE,CAL_STR
EXTERNAL::ELE_TRANS,FORM_K
OPEN(UNIT=1,FILE='D:\MY PROGRAM IN FORTRAN\A\BASIC.TXT')
OPEN(UNIT=8,FILE='D:\MY PROGRAM IN FORTRAN\A\DATA.TXT')
!OPEN(UNIT=10,FILE='D:\MY PROGRAM IN FORTRAN\A\TEST.TXT')
READ(1,*) ID,N_NODE,N_ELE,N_INT,PE,PT,PR,ROU
5 ,N_BC,N_LOAD,N_LINE_LOAD,X_N_LOAD,N_LINEAR(1:2),Q_LINEAR
ALLOCATE(G_ELE(1:N_ELE,8),MARK(1:N_ELE),G_NODE(1:N_NODE,2))
ALLOCATE(LINE_LOAD_NODE(1:N_LINE_LOAD,2),P_LOAD(1:N_LOAD,2))
ALLOCATE(X_LOAD_NODE(1:X_N_LOAD,2),X_LOAD(1:X_N_LOAD))
ALLOCATE(NP_LOAD(1:N_LOAD),LINE_LOAD(1:N_LINE_LOAD))
ALLOCATE(N_DIS(1:N_BC),BC(1:N_BC,2))
ALLOCATE(VV(1:N_NODE*2),DIAG(2*N_NODE+1),PP(2*N_NODE))
ALLOCATE(SIGMA(1:N_ELE,3))
ALLOCATE(AREA(1:N_ELE))
!READ DATA FROM DATA.TXT
CALL READ_DATA()
PRINT*,'DONE WITH READ_DATA'
CALL CAL_PE()
PRINT*,'DONE WITH CAL_PE'
!FORM THE MATRIX OF ELEMENT K
!CALL FORM_KE(S34,N_INT,PT,PR,PE,COORD,KE)
!FORM THE MATRIX OF TOTAL K
CALL BAND_K()
PRINT*,'DONE WITH BAND_K'
ALLOCATE(KK(DIAG(2*N_NODE+1)-1))
KK=0.0
DOF = DIAG(2*N_NODE+1)-1
CALL FORM_K()
PRINT*,'DONE WITH FORM_K'
!FORM THE SCALAR OF LOAD
VV=0.0
CALL FORM_LOAD()
PRINT*,'DONE WITH FORM_LOAD'
!ADD THE BOUNDARY CONDITION OF DISPLACEMENT
CALL ADD_BC()
PRINT*,'DONE WITH ADD_BC'
!SOLVE THE EQUATION
CALL SOLVE(KK,PP,DIAG,VV,N_NODE*2)
PRINT*,'DONE WITH SOLVE'
!CALCULATE THE STRESS
CALL CAL_STR()
PRINT*,'DONE WITH CAL_STR'
WRITE(8,*) 'DISPLACEMENT OF EACH NODE'
DO I = 1,2*N_NODE,2
WRITE(8,*) (I+1)/2,':',VV(I),VV(I+1)
END DO
WRITE(8,*)
WRITE(8,*) 'SIGMA FOR EACH NODE AFTER AVERAGING'
DO I = 1,N_ELE
WRITE(8,*) I,':',SIGMA(I,1),SIGMA(I,1),SIGMA(I,3)
END DO
WRITE(8,*)
PRINT*,'DONE WITH WRITING OUTPUT'
CLOSE(1)
CLOSE(6)
CLOSE(8)
END PROGRAM MAIN
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -