betterr_inter_sample.zip
源代码：
Module global
Implicit None
Integer, Save :: elem(1000000, 8), pindex(1000000), pmatche(100000), ii, elnum, nnum, pnum, kk, iii, key, ninnitial, keya, vec(3, 8), count, courier, e(1000000)
!单元节点,插值点编号索引,插值点与单元对应关系
Real (Kind=8), Save :: node(1000000, 3), erange(1000000, 6), point(100000, 3), locc(100000, 3), innitial(100, 3), toler, error(3), judge(2, 3), l(3), s(8), ds(8, 3), f(3), jacob(3, 3), dl(3), lamita, fsave(3, 2), lsave(3)
!节点坐标，单元的坐标范围，插值点坐标,L代表本地坐标,S(8)代表8个型函数，DS(8,3)代表八个型函数在不同自变量方向上的偏导数
End Module global
!*********************************主程序***********************************************
Program main
Use global
Implicit None
Call intro
Call input
Do ii = 1, pnum
iii = pindex(ii)
Do kk = 1, elnum
Call test
If (key==1) Then
Exit
End If
End Do
End Do
Pause
End Program main
!******************************输入数据的子程序************************************************************
Subroutine input
Use global
Integer i, k, j, a
Write (*, *) 'READING FILES'
Open (20, File='MESH.DAT')
Open (30, File='POINT.DAT')
Open (50, File='LOC.DAT')
Open (60, File='UFOUND.DAT')
Open (70, File='DEBUGINFOR.TXT')
Open (80, File='LOCCOD.TXT')
Open (170, File='NOTE.TXT')
Read (20, *) elnum, nnum, ninnitial
Do i = 1, nnum
Read (20, *) a, node(a, :)
End Do
Do i = 1, elnum
Read (20, *) a, elem(a, :)
e(i) = a
End Do
Do i = 1, ninnitial
Read (20, *) a, innitial(i, :)
End Do
Read (30, *) pnum, toler, error
Do i = 1, pnum
Read (30, *) pindex(i), (point(pindex(i),k), k=1, 3)
End Do
Close (20)
Close (30) !开始做erange数组
Write (50, *) pnum
Do i = 1, elnum
erange(e(i), 1) = minval(node(elem(e(i),:),1)) !要了解min（）函数与minval()函数的关系.min函数括号里面跟的是变量，minval跟的是数组。同时，隐循环不能出现在表达式中，只能出现在IO列表中。
erange(e(i), 2) = maxval(node(elem(e(i),:),1))
erange(e(i), 3) = minval(node(elem(e(i),:),2))
erange(e(i), 4) = maxval(node(elem(e(i),:),2))
erange(e(i), 5) = minval(node(elem(e(i),:),3))
erange(e(i), 6) = maxval(node(elem(e(i),:),3))
End Do
Write (*, *) 'FINISH READING'
End Subroutine input
!*********************寻找本地坐标子程序********************************
Subroutine test
Use global
key = 0
If (point(iii,1)>=erange(e(kk),1)-0.1 .And. point(iii,1)<=erange(e(kk),2)+0.1 .And. point(iii,2)>=erange(e(kk),3)-0.1 .And. point(iii,2)<=erange(e(kk),4)+0.1 .And. point(iii,3)>=erange(e(kk),5)-0.1 .And. point(iii,3)<=erange(e(kk),6)+0.1) Then
Call issac_newton
If (key==1) Then
Return
End If
If (e(kk)==elnum) Then
Write (*, *) 'POINT', iii, 'NOT FOUND'
Write (60, *) 'POINT', iii, 'NOT FOUND'
End If
End If
End Subroutine test
!*********************牛顿爵士的迭代法********************************
Subroutine issac_newton
Include 'link_f90_static.h'
Use wrrrn_int
Use global
Implicit None
Integer i, j, k
Real :: monitor1, monitor2
key = 0
count = 1
j = 0
l = innitial(count, :)
!JUDGE=20000.0
Do While (.True.)
j = j + 1
courier = 0
jacob = 0.0
lamita = 1
Call fem !做雅可比矩阵
fsave(:, 1) = f
lsave = l
monitor1 = abs(fsave(1,1)) + abs(fsave(2,1)) + abs(fsave(3,1))
Do While (.True.) !这个循环必须新值比旧值小才能跳出，符合牛顿下山的精神
Call fem
f = f
fsave(:, 1) = f
monitor1 = abs(fsave(1,1)) + abs(fsave(2,1)) + abs(fsave(3,1)) !MONITOR1是无穷向量
!CALL WRRRN('jacob',jacob)
Call solve(jacob, -1*f, dl, courier)
!CALL WRRRN('JACOB',JACOB)
l = l + dl
Call fem
fsave(:, 2) = f
monitor2 = abs(fsave(1,2)) + abs(fsave(2,2)) + abs(fsave(3,2))
If (monitor2<=0.1) Then !如果新F的无穷向量小的话，就直接跳出循环。
Exit
End If
If (monitor1>=monitor2) Then
Exit
End If
!LAMITA=LAMITA/2
!IF(LAMITA<=0.000001)THEN
! WRITE(170,*)'初值选在了坑里面'
!COUNT=COUNT+1
!L=INNITIAL(COUNT,:)
!J=0
!KEYA=0
!JUDGE=20000.0
!EXIT
!ENDIF
!L=LSAVE
End Do
!IF(J==0)THEN
! CYCLE
! ENDIF
!!!WRITE(*,'(a12,I8,a12,I8,a12,I8)')'插值点为',II,'单元为',e(KK),'初值更换次数为',COUNT
Write (170, *)
Write (170, *)
Write (170, '(a8,I8,A10,I8,A18,I8,A14,I8)') '插值点为', ii, '单元为', e(kk), '初值更换次数为', count, '迭代步数为', j
!!!WRITE(*,'(a12,I8)')'迭代步数为',J
!!!WRITE(*,'(a10,3F20.8)')'迭代增量为',DL
Write (170, '(a10,3F20.8)') '目标量为', l
Write (170, '(a10,3F20.8)') '迭代增量为', dl
!!!WRITE(*,*)
!IF((-1)**J==-1)THEN
! JUDGE(1,:)=DL
!ELSE IF((-1)**J==1)THEN
! JUDGE(2,:)=DL
!ENDIF
!
!
!IF ( ABS((JUDGE(1,1)+JUDGE(2,1)))<=0.000001 .AND. ABS((JUDGE(1,2)+JUDGE(2,2))<=0.000001) .AND. ABS((JUDGE(1,3)+JUDGE(2,3))<=0.000001))THEN
!KEYA=1
!ENDIF
If (courier==1 .Or. abs(dl(1))>=10000 .Or. abs(dl(2))>=10000 .Or. abs(dl(3))>=10000 .Or. j>=1000 .Or. keya==1) Then
!给出了判断不收敛的4个条件：1.迭代方程无唯一解 2.解的数值过大 3.迭代次数过多 4.发生数值震荡。 这4个条件是在调程序时候遇到的问题的解决办法。
If (count==ninnitial) Then
Write (170, *) '穷尽初值，未能收敛，请增加初值。'
Write (170, *) '****************************************************************'
Return
End If
!!!WRITE(*,*)'迭代不收敛，更换初值再次迭代'
Write (170, *) '迭代不收敛，更换初值再次迭代'
Write (170, *) '****************************************************************'
count = count + 1
l = innitial(count, :)
j = 0
keya = 0
!JUDGE=20000.0
Cycle
End If
If (monitor2<=0.01) Then
!!!write(*,'(a12,I8,A12,3F20.8)')'收敛迭代次数',j,'本地坐标为',l
Write (170, '(a12,I8,A12,3F20.8)') '收敛迭代次数', j, '本地坐标为', l
Write (170, *) '****************************************************'
Write (170, *)
Write (170, *)
Write (70, '(4I8,3F20.8)') iii, e(kk), j, l
If (abs(l(1))<=1.00+error(1) .And. abs(l(2))<=1.00+error(2) .And. abs(l(3))<=1.00+error(3)) Then
locc(iii, :) = l
Write (80, '(2i8,3f20.8)') ii, e(kk), l
Write (*, *) '第', ii, '个点的本地坐标已经找到了'
key = 1
End If
Return
End If
End Do
End Subroutine issac_newton
!*********************高斯消元法********************************
Subroutine solve(left, right, x, courier)
Use wrrrn_int
Integer i, j, courier
Real (Kind=8) :: left(3, 3), right(3), a(3, 4), b(3), c, d(4), x(3)
a(1, 1:3) = left(1, :)
a(2, 1:3) = left(2, :)
a(3, 1:3) = left(3, :)
a(:, 4) = right
c = maxval(abs(a(:,1)))
Do i = 1, 3
If (abs(a(i,1))==c) Then
Goto 100
End If
End Do
100 Continue
d = a(1, :)
a(1, :) = a(i, :)
a(i, :) = d
!CALL WRRRN('A',A)
Do i = 2, 3
c = a(i, 1)
Do j = 1, 4
a(i, j) = a(1, j)*c - a(i, j)*a(1, 1)
End Do
End Do
!CALL WRRRN('A',A)
c = maxval(abs(a(2:3,2)))
Do i = 2, 3
If (abs(a(i,2))==c) Then
Goto 300
End If
End Do
300 Continue
d = a(2, :)
a(2, :) = a(i, :)
a(i, :) = d
c = a(3, 2)
Do j = 2, 4
a(3, j) = a(2, j)*c - a(3, j)*a(2, 2)
End Do
!CALL WRRRN('A',A)
If (a(1,1)*a(2,2)*a(3,3)==0) Then
courier = 1
Return
End If
x(3) = a(3, 4)/a(3, 3)
x(2) = (a(2,4)-a(2,3)*x(3))/a(2, 2)
x(1) = (a(1,4)-a(1,3)*x(3)-a(1,2)*x(2))/a(1, 1)
End Subroutine solve
!*********************高斯消元法********************************
Subroutine fem
Use global
Integer i, j, k
f = 0.0
Do i = 1, 8
vec(1, i) = (-1)**floor(1.5+0.45*i) !系数
vec(2, i) = (-1)**floor(1.05+0.45*i)
vec(3, i) = (-1)**floor(1.05+0.20*i)
s(i) = 0.125*(1+vec(1,i)*l(1))*(1+vec(2,i)*l(2))*(1+vec(3,i)*l(3)) !型函数N（i）
ds(i, 1) = 0.125*vec(1, i)*(1+vec(2,i)*l(2))*(1+vec(3,i)*l(3)) !型函数n(i)的偏导
ds(i, 2) = 0.125*vec(2, i)*(1+vec(1,i)*l(1))*(1+vec(3,i)*l(3))
ds(i, 3) = 0.125*vec(3, i)*(1+vec(1,i)*l(1))*(1+vec(2,i)*l(2))
f(1) = s(i)*node(elem(e(kk),i), 1) + f(1)
f(2) = s(i)*node(elem(e(kk),i), 2) + f(2)
f(3) = s(i)*node(elem(e(kk),i), 3) + f(3)
If (i==8) Then
f(1) = f(1) - point(iii, 1)
f(2) = f(2) - point(iii, 2)
f(3) = f(3) - point(iii, 3)
End If
jacob(1, 1) = ds(i, 1)*node(elem(e(kk),i), 1) + jacob(1, 1) !雅可比矩阵一定要仔细想想，做对
jacob(1, 2) = ds(i, 2)*node(elem(e(kk),i), 1) + jacob(1, 2)
jacob(1, 3) = ds(i, 3)*node(elem(e(kk),i), 1) + jacob(1, 3)
jacob(2, 1) = ds(i, 1)*node(elem(e(kk),i), 2) + jacob(2, 1)
jacob(2, 2) = ds(i, 2)*node(elem(e(kk),i), 2) + jacob(2, 2)
jacob(2, 3) = ds(i, 3)*node(elem(e(kk),i), 2) + jacob(2, 3)
jacob(3, 1) = ds(i, 1)*node(elem(e(kk),i), 3) + jacob(3, 1)
jacob(3, 2) = ds(i, 2)*node(elem(e(kk),i), 3) + jacob(3, 2)
jacob(3, 3) = ds(i, 3)*node(elem(e(kk),i), 3) + jacob(3, 3)
End Do
End Subroutine fem
Subroutine intro
Write (*, *) '*************************************************************'
Write (*, *) '* by：石中岳 *'
Write (*, *) '* tips：仅供学习使用，提供源代码 *'
Write (*, *) '* 欢迎讨论，联系方式qq641240961 *'
Write (*, *) '*************************************************************'
End Subroutine intro
