首页 >

有限元插值搜寻本地坐标程序

作者:石中岳  日期:03-20
来源:Fcode研讨团队
范例数据:
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

常规|工具|专业|读物|
代码|教学|算法|
首页 >
FortranCoder手机版-导航