📄 main.f90
字号:
!***************************************************************************!
! 主程序功能： !
! !
! 求解第一类边界条件、三角元剖分、线性插值的位场延拓 !
! 作者：马见青 !
! 时间：2008年7月 !
!---------------------------------------------------------------------------!
! 摘自 《地球物理中的有限单元法》(徐世浙 著) !
!***************************************************************************!
program FEM
character*40::filename1='i3_input.txt'
character*40::filename2='xy_input.txt'
character*40::filename3='nb1_input.txt'
character*40::filename4='u1_input.txt'
integer,parameter::nd=15 ! 节点总数
integer,parameter::ne=16 ! 单元总数
integer,parameter::nd1=12 ! 第一类边界节点数
integer i3(3,ne) ! 单元节点编号
real xy(2,nd) ! 节点坐标
integer nb1(nd1) ! 第一类边界节点号
real u1(nd1) ! 场值
real ke(3,3)
real,allocatable::sk(:,:)
real,allocatable::u(:)
!-------------从文件i3_input.txt中读取单元节点编号-----------------
write(*,*) '*******************************************************'
write(*,*) '单元节点编号:'
open(10,file=filename1)
do i=1,NE
read(10,*) I3(1,i),I3(2,i),I3(3,i)
write(*,*) I3(1,i),I3(2,i),I3(3,i)
end do
close(10)
!-------------从文件xy_input.txt中读取节点坐标---------------------
write(*,*) '*******************************************************'
write(*,*) '节点坐标:'
open(20,file='xy_input.txt')
do i=1,ND
read(20,*) XY(1,i),XY(2,i)
write(*,*) XY(1,i),XY(2,i)
end do
close(20)
!------------从文件nb1_input.txt中读取第一类边界节点号-------------
write(*,*) '*******************************************************'
write(*,*) '已知节点号:'
open(30,file='nb1_input.txt')
do i=1,nd1
read(30,*) nb1(i)
write(*,*) NB1(i)
end do
close(30)
!------------------从文件u1_input.txt中读取场值---------------------
write(*,*) '*******************************************************'
write(*,*) '已知节点号的场值:'
open(40,file='u1_input.txt')
do i=1,nd1
read(40,*) u1(i)
write(*,*) u1(i)
end do
close(40)
!===================================================================
write(*,*) '*******************************************************'
call MBW(ne,i3,iw)
write(*,*) '半带宽:',iw
!--------------------------------------------------------------------
write(*,*) '*******************************************************'
allocate(SK(ND,IW))
call UK1(nd,ne,iw,i3,xy,sk)
! call UK1(ND,NE,IW,I3,XY,SK)
write(*,*) '总体系数矩阵:'
do i=1,ND
write(*,*) (SK(i,j),j=1,IW)
end do
!--------------------------------------------------------------------
write(*,*) '*******************************************************'
allocate(u(nd))
call UB1(nd1,nb1,u1,nd,iw,sk,u)
write(*,*) '加入第一类边界条件后的总体系数矩阵:'
do i=1,ND
write(*,*) (SK(i,j),j=1,IW)
end do
write(*,*) '*******************************************************'
write(*,*) '加入第一类边界条件后的右端向量:'
do i=1,ND
write(*,*) U(i)
end do
!--------------------------------------------------------------------
write(*,*) '*******************************************************'
call LDLT(SK,ND,IW,U,IE)
write(*,*) '解向量:'
do i=1,ND
write(*,*) U(i)
end do
end program
!***************************************************************************!
! subroutine UKE1(x,y,ke) !
! 功能 三角单元、线性插值时计算单元系数矩阵ke。 !
! 参数说明 !
! x----3个元的一维实数组,输入参数,存放单元节点的x坐标,存放顺序xi,xj,xm. !
! y----3个元的一维实数组,输入参数,存放单元节点的y坐标,存放顺序yi,yj,ym. !
! ke---3*3的二维实数组,输出参数,存放单元系数矩阵. !
!***************************************************************************!
SUBROUTINE UKE1(X,Y,KE)
IMPLICIT NONE
REAL X(3)
REAL Y(3)
REAL KE(3,3)
REAL A(3)
REAL B(3)
INTEGER I,J
REAL S
A(1)=Y(2)-Y(3)
A(2)=Y(3)-Y(1)
A(3)=Y(1)-Y(2)
B(1)=X(3)-X(2)
B(2)=X(1)-X(3)
B(3)=X(2)-X(1)
S=2.0*(A(1)*B(2)-A(2)*B(1))
DO 10 I=1,3
DO 10 J=1,I
10 KE(I,J)=(A(I)*A(J)+B(I)*B(J))/S
RETURN
END SUBROUTINE UKE1
!***************************************************************************!
! SUBROUTINE MBW(ne,i3,iw) !
! !
! 功能 计算总体系数矩阵的半带宽。 !
! 参数说明 !
! ne---整型变量,输入参数,单元总数。 !
! i3---3*ne的二维整数组,输入参数,存放单元节点编号。 !
! iw---整型变量,输出参数,半带宽。 !
!***************************************************************************!
SUBROUTINE MBW(NE,I3,IW)
IMPLICIT NONE
INTEGER NE
INTEGER I3(3,NE)
INTEGER IW
INTEGER I
INTEGER M
DO I=1,NE
M=MAX(IABS(I3(1,I)-I3(2,I)),IABS(I3(2,I)-I3(3,I)),IABS(I3(3,I)-I3(1,I)))
IF(M+1>IW) IW=M+1
END DO
RETURN
END SUBROUTINE MBW
!***************************************************************************!
! SUBROUTINE UK1(nd,ne,iw,i3,xy,sk) !
! !
! 功能 三角单元,线性插值时用定带宽储存的方法集成总体矩阵 !
! 参数说明 !
! nd---整型变量,输入参数,节点总数。 !
! ne---整型变量,输入参数,单元总数。 !
! iw---整型变量,输入参数,半带宽。 !
! i3---3*ne的二维整数组,输入参数,存放单元节点编号,存放顺序:i1,j1,m1, !
! ...,ine,jne,mne。 !
! xy---2*nd的二维整数组,输入参数,存放节点的xy坐标,存放顺序:x1,y1, !
! ...,xnd,ynd。 !
! sk---nd*iw的二维整数组,输出参数,存放定带宽储存的总体系数矩阵 !
!***************************************************************************!
SUBROUTINE UK1(ND,NE,IW,I3,XY,SK)
IMPLICIT NONE
INTEGER ND
INTEGER NE
INTEGER IW
INTEGER I3(3,NE)
REAL XY(2,ND)
REAL SK(ND,IW)
REAL X(3)
REAL Y(3)
REAL KE(3,3)
INTEGER I,J,K,L
INTEGER NJ
INTEGER NK
DO 10 I=1,ND
DO 10 J=1,IW
10 SK(I,J)=0.
DO 20 L=1,NE
DO 30 J=1,3
I=I3(J,L)
X(J)=XY(1,I)
30 Y(J)=XY(2,I)
CALL UKE1(X,Y,KE)
DO 40 J=1,3
NJ=I3(J,L)
DO 40 K=1,J
NK=I3(K,L)
IF(NJ<NK) GOTO 50
NK=NK-NJ+IW
SK(NJ,NK)=SK(NJ,NK)+KE(J,K)
GOTO 40
50 NJ=NJ-NK+IW
SK(NK,NJ)=SK(NK,NJ)+KE(J,K)
NJ=NJ+NK-IW
40 CONTINUE
20 CONTINUE
RETURN
END SUBROUTINE UK1
!***************************************************************************!
! SUBROUTINE UB1(nd1,nb1,u1,nd,iw,sk,u) !
! 功能 在定带宽储存的总体系数矩阵和右侧列向量上加上第一类边界条件。 !
! 参数说明 !
! nd1---整型变量,输入参数,第一类边界条件的节点数。 !
! nb1---nd1个元的一维整数组，输入参数，存放第一类边界条件的节点号。 !
! u1----nd1个元的一维实数组，输入参数，存放第一类边界节点的场值。 !
! nd----整型变量,输入参数,节点总数。 !
! iw----整型变量,输入参数,半带宽。 !
! sk----nd*iw的二维整数组,输入、输出参数,定带宽存放总体系数矩阵，输出 !
! 时,第一类边界条件已带入。 !
! u-----nd个元的一维实数组，输出参数，存放加入第一类边界条件后的右端列 !
! 向量。 !
!***************************************************************************!
SUBROUTINE UB1(ND1,NB1,U1,ND,IW,SK,U)
IMPLICIT NONE
INTEGER ND1
INTEGER ND
INTEGER IW
INTEGER NB1(ND1)
REAL U1(ND1)
REAL SK(ND,IW)
REAL U(ND)
INTEGER I,J
DO 10 I=1,ND
10 U(I)=0.0
DO 20 I=1,ND1
J=NB1(I)
SK(J,IW)=SK(J,IW)*1.0E10
20 U(J)=SK(J,IW)*U1(I)
RETURN
END SUBROUTINE UB1
!***************************************************************************!
! SUBROUTINE LDLT(a,n,iw,p,ie) !
! 功能 对称带型线性方程组的系数矩阵的下三角部分被定带宽地储存在矩型 !
! 数组a中，利用a数组来解方程组。 !
! 参数说明 !
! a-----n*iw的二维实数组,输入参数,存放对称带型线性方程组的系数矩阵的 !
! 下三角部分。 !
! n-----整型变量,输入参数,方程组的阶数。 !
! iw----整型变量,输入参数,对称带型线性方程组的半带宽。 !
! p-----n个元素的一维实数组，输入、输出参数。开始存放方程组的右侧列 !
! 向量，工作结束时，存放解向量。 !
! ie----整型变量，输出参数，标志。ie=0时表示程序正常结束； !
! ie=1时，表示系数矩阵奇异。 !
!***************************************************************************!
SUBROUTINE LDLT(A,N,IW,P,IE)
IMPLICIT NONE
INTEGER N
INTEGER IW
INTEGER IE
REAL A(N,IW)
REAL P(N)
INTEGER I,J,K,L,IT,MI,MJ
INTEGER IL,IJ,JL
REAL B
DO 15 I=1,N
IF(I<=IW) GOTO 20
IT=I-IW+1
GOTO 30
20 IT=1
30 K=I-1
IF(I==1) GOTO 40
DO 25 L=IT,K
IL=L+IW-I
B=A(I,IL)
A(I,IL)=B/A(L,IW)
P(I)=P(I)-A(I,IL)*P(L)
MI=L+1
DO 25 J=MI,I
IJ=J+IW-I
JL=L+IW-J
25 A(I,IJ)=A(I,IJ)-A(J,JL)*B
40 IF(A(I,IW)==0.0) GOTO 100 !/ ??????????????
15 CONTINUE
DO 45 J=1,N
IF(J<=IW) GOTO 60
IT=N-J+IW
GOTO 70
60 IT=N
70 I=N-J+1
P(I)=P(I)/A(I,IW)
IF(J==1) GOTO 45
K=I+1
DO 65 MJ=K,IT
IJ=I-MJ+IW
65 P(I)=P(I)-P(MJ)*A(MJ,IJ)
45 CONTINUE
IE=0
GOTO 110
100 IE=1
110 RETURN
END SUBROUTINE LDLT
!--------------------------------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -