C Solution of -Δu = f with Dirichlet Boundary Condition
C by Finite Element Method on the rectangle [0,1]x[0,1]
C 章末課題 12.3 に対応するためグラフィック表示をやめ, 分割数を
C 自動的に 10, 20, 40, 80 と取って Richardson 補外法により真値を
C 予測し, 誤差のオーダーを見るように変更したもの
C 値の比較は中央値だけで行っている
C
C    640 25281 51200 Band Breadth : 160
C   1280    101761    204800
      PROGRAM POISN3
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000,IBWMAX=1000)
      PARAMETER(IRCT=5)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO),
     +       A(0:IBWMAX,NIPNO),U(NIPNO),V(NIPNO),T(0:IRCT-1,IRCT)
      INTEGER*4 ITRI(ITRINO)

C      WRITE(*,*)'Division Number of X, Y axis Resp. ( =< 100) : '
C      READ(*,*)NX,NY
C
      X0=0.0D0
      X1=1.0D0
      Y0=0.0D0
      Y1=1.0D0
      NX=10
      NY=10
      DO 10000 KAISU=1,IRCT
C Manual Generation of Triangulation
         CALL TRIANG(VERTX,VERTY,ITRI,NBP,NIP,NT,
     +            X0,X1,Y0,Y1,NX,NY,HX,HY)
C Preparation of Finite Approximatory Matrices
         CALL PREPMV(A,V,VERTX,VERTY,ITRI,NIP,NT,IBW)
         WRITE(*,'(1H ,A14,I4)')'Band Breadth :',IBW
C Resolution of Approximate Linear Equation
         CALL SOLVLE(A,U,V,NIP,IBW)
C Drawing the Graph of the Result
C      CALL CLS()
C      CALL DRWSOL(U,VERTX,VERTY,ITRI,NT,X0,X1,Y0,Y1,
C     +            -1.0D0,1.0D0)
C Value at the ASKING THE VALUE OF THE APPROXIMATE SOLUTION
C         Z=VALUE(0.5D0,0.5D0,U,VERTX,VERTY,ITRI,NT,HX,HY)
         Z=U(NIP/2+1)
         WRITE(*,*)'Value at the center : ',Z
         T(0,KAISU)=Z
         NX=NX*2
         NY=NY*2
10000 CONTINUE
      M=2
      DO 300 I=1,IRCT-1
         DO 200 I2=1,IRCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
  200    CONTINUE
         M=M*2
  300 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 1st order:'
      DO 500 I=0,IRCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
  500 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IRCT-1,1)),I=1,IRCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IRCT-1,1))/
     +            (T(0,I)-T(IRCT-1,1)),I=1,IRCT-1)
      M=4
      DO 3000 I=1,IRCT-1
         DO 2000 I2=1,IRCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 2000    CONTINUE
         M=M*4
 3000 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 2nd order:'
      DO 5000 I=0,IRCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 5000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IRCT-1,1)),I=1,IRCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IRCT-1,1))/
     +              (T(0,I)-T(IRCT-1,1)),I=1,IRCT-1)

      END
C
      SUBROUTINE TRIANG(VERTX,VERTY,ITRI,NBP,NIP,NT,
     +                  X0,X1,Y0,Y1,NX,NY,HX,HY)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000)
C      PARAMETER(NBPNO=-1000,NIPNO=30000,ITRINO=100000)
C      PARAMETER(NBPNO=-800,NIPNO=10000,ITRINO=60000)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO)
      INTEGER*4 ITRI(ITRINO),ITRIS(ITRINO/30)
      CHARACTER CH
C 初期化
      NBP=0
      NIP=0
      NT=0
      IBW=0
      NT3=0
      HX=(X1-X0)/NX
      HY=(Y1-Y0)/NY
      X=X0
      Y=Y0
C 左下角の第１三角形の登録
      NT=NT+1
      NBP=NBP-1
      VERTX(NBP)=X
      VERTY(NBP)=Y
      NT3=NT3+1
      ITRI(NT3)=NBP
      NIP=NIP+1
      VERTX(NIP)=X+HX
      VERTY(NIP)=Y+HY
      NT3=NT3+1
      ITRI(NT3)=NIP
      NBP=NBP-1
      VERTX(NBP)=X
      VERTY(NBP)=Y+HY
      NT3=NT3+1
      ITRI(NT3)=NBP
C 左下角の第２三角形の登録
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(1)
      NBP=NBP-1
      VERTX(NBP)=X+HX
      VERTY(NBP)=Y
      NT3=NT3+1
      ITRI(NT3)=NBP
      NT3=NT3+1
      ITRI(NT3)=ITRI(2)
C 最下行の三角形の登録
      DO 100 I=2,NX-1
          X=X+HX
          NT=NT+1
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-2)
          NT3=NT3+1
          NIP=NIP+1
          VERTX(NIP)=X+HX
          VERTY(NIP)=Y+HY
          ITRI(NT3)=NIP
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-3)
          NT=NT+1
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-3)
          NT3=NT3+1
          NBP=NBP-1
          VERTX(NBP)=X+HX
          VERTY(NBP)=Y
          ITRI(NT3)=NBP
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-4)
  100 CONTINUE
C 右下角の第１三角形の登録
      X=X+HX
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-2)
      NBP=NBP-1
      VERTX(NBP)=X+HX
      VERTY(NBP)=Y
      NT3=NT3+1
      ITRI(NT3)=NBP
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-3)
C 右下角の第２三角形の登録
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-2)
      NBP=NBP-1
      VERTX(NBP)=X+HX
      VERTY(NBP)=Y+HY
      NT3=NT3+1
      ITRI(NT3)=NBP
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-3)
C 中央の三角形の登録
      DO 250 J=2,NY-1
         Y=Y+HY
         X=X0
C 左の第１三角形の登録
         NT=NT+1
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-6*NX+2)
         NT3=NT3+1
         NIP=NIP+1
         VERTX(NIP)=X+HX
         VERTY(NIP)=Y+HY
         ITRI(NT3)=NIP
         NT3=NT3+1
         NBP=NBP-1
         VERTX(NBP)=X
         VERTY(NBP)=Y+HX
         ITRI(NT3)=NBP
C 左の第２三角形の登録
         NT=NT+1
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-3)
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-6*NX+1)
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-4)
C 中央の登録
         DO 200 I=2,NX-1
            X=X+HX
            NT=NT+1
            NT3=NT3+1
            ITRI(NT3)=ITRI(NT3-2)
            NT3=NT3+1
            NIP=NIP+1
            VERTX(NIP)=X+HX
            VERTY(NIP)=Y+HY
            ITRI(NT3)=NIP
            NT3=NT3+1
            ITRI(NT3)=ITRI(NT3-3)
            NT=NT+1
            NT3=NT3+1
            ITRI(NT3)=ITRI(NT3-3)
            NT3=NT3+1
            ITRI(NT3)=ITRI(NT3-6*NX+1)
            NT3=NT3+1
            ITRI(NT3)=ITRI(NT3-4)
  200    CONTINUE
C 右の第１三角形の登録
         X=X+HX
         NT=NT+1
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-2)
         NBP=NBP-1
         VERTX(NBP)=X+HX
         VERTY(NBP)=Y+HY
         NT3=NT3+1
         ITRI(NT3)=NBP
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-3)
C 右の第２三角形の登録
         NT=NT+1
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-3)
         NT3=NT3+1
         IF (J.GT.2) THEN
            ITRI(NT3)=ITRI(NT3-6*NX+1)
         ELSE
            ITRI(NT3)=ITRI(NT3-6*NX)
         END IF
         NT3=NT3+1
         ITRI(NT3)=ITRI(NT3-4)
  250 CONTINUE    
C 左上角の第１三角形の登録
      Y=Y+HY
      X=X0
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-6*NX+2)
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-6*NX)
      NT3=NT3+1
      NBP=NBP-1
      VERTX(NBP)=X
      VERTY(NBP)=Y+HY
      ITRI(NT3)=NBP
C 左上角の第２三角形の登録
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-2)
      NBP=NBP-1
      VERTX(NBP)=X+HX
      VERTY(NBP)=Y+HY
      NT3=NT3+1
      ITRI(NT3)=NBP
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-3)
C 左上の第３三角形の登録
      X=X+HX
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-3)
      NT3=NT3+1
      NBP=NBP-1
      VERTX(NBP)=X+HX
      VERTY(NBP)=Y+HY
      ITRI(NT3)=NBP
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-4)
C 左上の第４三角形の登録
      NT=NT+1
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-3)
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-6*NX+1)
      NT3=NT3+1
      ITRI(NT3)=ITRI(NT3-4)
C 最上行の三角形の登録
      DO 300 I=3,NX
          X=X+HX
          NT=NT+1
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-2)
          NT3=NT3+1
          NBP=NBP-1
          VERTX(NBP)=X+HX
          VERTY(NBP)=Y+HY
          ITRI(NT3)=NBP
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-3)
          NT=NT+1
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-3)
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-6*NX+1)
          NT3=NT3+1
          ITRI(NT3)=ITRI(NT3-4)
  300 CONTINUE
C      CALL DRWTRI(VERTX,VERTY,ITRI,NT,X0,X1,Y0,Y1)
      WRITE(*,*)'NUMBER OF '
      WRITE(*,*)'BOUNDARY VERTICES, INNER VERTICES, TRIANGLES :'
      WRITE(*,'(1H ,3I10)')-NBP,NIP,NT
C 3600 WRITE(*,*)'CONTINUE? Y/N : '
C      READ(*,'(A1)')CH
C      IF (INDEX('Yy',CH).EQ.0) STOP
      RETURN
      END
C
C Draws the Result of the Triangulation
      SUBROUTINE DRWTRI(VERTX,VERTY,ITRI,NT,X0,X1,Y0,Y1)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000)
C      PARAMETER(NBPNO=-1000,NIPNO=30000,ITRINO=100000)
C      PARAMETER(NBPNO=-800,NIPNO=10000,ITRINO=60000)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO)
      INTEGER*4 ITRI(ITRINO),II(3),IX(3),IY(3)
      PARAMETER(IXMIN=0,IXMAX=799,IYMIN=0,IYMAX=599,IR=8)
      IGX(X)=IXMAX-IYMAX+(X-XMIN)/(XMAX-XMIN)*(IYMAX-IYMIN)
      IGY(Y)=IYMIN+(Y-YMIN)/(YMAX-YMIN)*(IYMAX-IYMIN)
      XH=(X1-X0)/20
      YH=(Y1-Y0)/20
      IF (XH.GE.YH) THEN
         XMIN=X0-XH
         XMAX=X1+XH
         YMIN=Y0-XH
         YMAX=Y0+X1-X0+XH
      ELSE
         YMIN=Y0-YH
         YMAX=Y1+YH
         XMAX=X1+YH
         XMIN=X1-Y1+Y0-YH
      END IF
C      CALL INIT(IXMIN,IYMIN,IXMAX,IYMAX)
C      CALL CLS()
      DO 1000 I=1,NT
        NT3=(I-1)*3
        DO 100 J=1,3
            II(J)=ITRI(NT3+J)
            IX(J)=IGX(VERTX(II(J)))
            IY(J)=IGY(VERTY(II(J)))
  100   CONTINUE
        DO 150 J=1,3
            IF (II(J).LT.0) THEN
C               CALL BLCOLR(5)
C               CALL CIRCLE(IX(J),IY(J),IR)
C            ELSE
C               CALL BLCOLR(9)
C               CALL BOX(IX(J)-IR,IY(J)-IR,IX(J)+IR,IY(J)+IR)
            END IF
  150   CONTINUE
        DO 200 J=1,3
            JJ=MOD(J,3)+1
C            CALL LINE(IX(J),IY(J),IX(JJ),IY(JJ))
  200   CONTINUE
 1000 CONTINUE
      RETURN
      END
C
C Preparation of the Stiffness Matrix and the Weight Vector.
C IBW retains the  breadth of the band matrix.
      SUBROUTINE PREPMV(A,V,VERTX,VERTY,ITRI,NIP,NT,IBW)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000,IBWMAX=1000)
C      PARAMETER(NBPNO=-1000,NIPNO=30000,ITRINO=100000,IBWMAX=200)
C      PARAMETER(NBPNO=-800,NIPNO=10000,ITRINO=60000,IBWMAX=100)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO),
     +       A(0:IBWMAX,NIPNO),V(NIPNO),X(3),Y(3),
     +       XX(3),YY(3)
      INTEGER*4 ITRI(ITRINO),II(3),IN(3)
      DO 100 J=1,NIP
        DO 90 I=0,IBWMAX
          A(I,J)=0
   90   CONTINUE
      V(J)=0
  100 CONTINUE
      IBW=0
      DO 1000 I=1,NT
       NT3=(I-1)*3
       DO 200 J=1,3
         II(J)=ITRI(NT3+J)
         IF (II(J).GT.0) IN(J)=II(J)
         X(J)=VERTX(II(J))
         Y(J)=VERTY(II(J))
  200  CONTINUE
       D=ORIENT(X(1),Y(1),X(2),Y(2),X(3),Y(3))
       DO 300 J=1,3
         J1=MOD(J,3)+1
         J2=MOD(J1,3)+1
         XX(J)=(X(J1)+X(J2))/2
         YY(J)=(Y(J1)+Y(J2))/2
  300  CONTINUE
       DO 400 J=1,3
         J1=MOD(J,3)+1
         J2=MOD(J1,3)+1
         IF (II(J).GT.0) THEN
            A(0,IN(J))=A(0,IN(J))
     +                     +((X(J1)-X(J2))*(X(J1)-X(J2))
     +                       +(Y(J1)-Y(J2))*(Y(J1)-Y(J2)))/2/D
            V(IN(J))=V(IN(J))+(F(XX(J1),YY(J1))+F(XX(J2),YY(J2)))*D/12
         END IF
         IF ((II(J1).GT.0).AND.(II(J2).GT.0)) THEN
            IDF=IABS(IN(J1)-IN(J2))
            IBW=MAX(IBW,IDF)
            OFFD=((X(J)-X(J1))*(X(J)-X(J2))
     +           +(Y(J)-Y(J1))*(Y(J)-Y(J2)))/2/D
            ISC=MIN(IN(J1),IN(J2))
            A(IDF,ISC)=A(IDF,ISC)-OFFD
         END IF
  400  CONTINUE
 1000 CONTINUE
      RETURN
      END
C
C Inner product of two vectors
C (Subroutine for Conjugate Gradient Method)
      FUNCTION AINPRD(V,W,N)
      REAL*8 AINPRD,V(N),W(N)
      AINPRD=0.0D0
      DO 100 I=1,N
          AINPRD=AINPRD+V(I)*W(I)
  100 CONTINUE
      RETURN
      END
C                     -1
C Calculation of U = A  V by Conjugate Gradient Method
      SUBROUTINE SOLVLE(A,U,V,N,IB)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(IBWMAX=1000,NIPNO=500000,EPSLON=1.0D-7)
C      PARAMETER(IBWMAX=100,NIPNO=10000,EPSLON=1.0D-7)
C THE CONSTANT EPSLON IS USED TO JUDGE THE CONVERGENCE
      REAL*8 A(0:IBWMAX,NIPNO),U(NIPNO),V(NIPNO),W0(NIPNO),W1(NIPNO),
     +       WR(NIPNO)
      DO 5 I=1,N
         WR(I)=V(I)
         W0(I)=V(I)
    5 CONTINUE
      RIPC=AINPRD(WR,WR,N)
      DO 500 KCONV=1,N-1
         DO 130 I=1,N
            W1(I)=0
            DO 110 J=MAX(1,I-IB),I
               W1(I)=W1(I)+A(I-J,J)*W0(J)
  110       CONTINUE
            DO 120 J=I+1,MIN(N,I+IB)
               W1(I)=W1(I)+A(J-I,I)*W0(J)
  120       CONTINUE
  130    CONTINUE
         AVIP=AINPRD(W0,W1,N)
         ALPHA=RIPC/AVIP
         DO 200 I=1,N
            U(I)=U(I)+ALPHA*W0(I)
            WR(I)=WR(I)-ALPHA*W1(I)
  200    CONTINUE
         RIP=AINPRD(WR,WR,N)
C Convergence test
         IF (RIP.LE.EPSLON) GO TO 1000
         BETA=RIP/RIPC
         RIPC=RIP
         DO 300 I=1,N
            W0(I)=WR(I)+BETA*W0(I)
  300    CONTINUE
  500 CONTINUE
 1000 RETURN
      END
C
C Draws the Graph of the Solution as a Surface in R^3
      SUBROUTINE DRWSOL(U,VERTX,VERTY,ITRI,NT,
     +                  X0,X1,Y0,Y1,Z0,Z1)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(IXMIN=0,IXMAX=799,IYMIN=0,IYMAX=599,IXD=5,IYD=6)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000)
C      PARAMETER(NBPNO=-1000,NIPNO=30000,ITRINO=100000)
C      PARAMETER(NBPNO=-800,NIPNO=10000,ITRINO=60000)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO),U(NIPNO),
     +       XW(3),YW(3),ZW(3)
      INTEGER*4 ITRI(ITRINO),II(3)
      CHARACTER*2 CHARX,CHARY,CHARZ
      IGX(X,Y,Z)=IXBIAS+IBREAD*(Y-Y0)/(Y1-Y0)
     +                 +ISLANT*(X1-X)/(X1-X0)
      IGY(X,Y,Z)=IYBIAS+ISLANT*(X1-X)/(X1-X0)
     +           +IBREAD*(Z-Z0)/(Z1-Z0)
C      IGX(X,Y,Z)=IXBIAS+AINCX*((Y-Y0)/(Y1-Y0)+(X1-X)/(X1-X0))
C      IGY(X,Y,Z)=INCY*((X1-X)/(X1-X0)
C     +           +(Y1-Y)/(Y1-Y0))+INCZ*(Z-Z0)/(Z1-Z0)
      PERSP=0.8
      IBREAD=(IYMAX-IYMIN+1)/(1.0+PERSP/1.414)
      ISLANT=IBREAD*PERSP/1.414
      IXBIAS=(IXMAX-IXMIN+1)-IBREAD-ISLANT-116
      IYBIAS=-116
C      IXBIAS=IXMAX/3
C      INCZ=(IYMAX+1)/2
C      AINCX=SQRT(3.0)*INCZ/2
C      INCY=INCZ/2
      CHARX='x'//CHAR(0)
      CHARY='y'//CHAR(0)
      CHARZ='z'//CHAR(0)
      XL=(X1-X0)/4
      YL=(Y1-Y0)/4
      ZL=(Z1-Z0)/4
C      CALL LCOLOR(11)
C      CALL CHCOLR(11)
      IX1=IGX(X0-XL,0.0D0,0.0D0)
      IY1=IGY(X0-XL,0.0D0,0.0D0)
      IX2=IGX(X1+XL,0.0D0,0.0D0)
      IY2=IGY(X1+XL,0.0D0,0.0D0)
C      CALL LINE(IX1,IY1,IX2,IY2)
      IXX=(IYMAX-IYMIN)/150
      IXX1=IXX*12/5
      IXX2=IXX*2
      IXX3=IXX*3
      IXX4=IXX*4
      IXX6=IXX*6
C      CALL LINE(IX2,IY2,IX2+IXX,IY2+IXX3)
C      CALL LINE(IX2,IY2,IX2+IXX3,IY2)
C      CALL CHARAC(IX2-IXX6,IY2-IXX2,CHARX)
      IX1=IGX(0.0D0,Y0-YL,0.0D0)
      IY1=IGY(0.0D0,Y0-YL,0.0D0)
      IX2=IGX(0.0D0,Y1+YL,0.0D0)
      IY2=IGY(0.0D0,Y1+YL,0.0D0)
C      CALL LINE(IX1,IY1,IX2,IY2)
C      CALL CHARAC(IX2+IXX2,IY2-IXX1,CHARY)
C      CALL LINE(IX2,IY2,IX2-IXX,IY2+IXX3)
C      CALL LINE(IX2,IY2,IX2-IXX3,IY2)
      IX1=IGX(0.0D0,0.0D0,Z0-ZL)
      IY1=IGY(0.0D0,0.0D0,Z0-ZL)
      IX2=IGX(0.0D0,0.0D0,Z1+ZL)
      IY2=IGY(0.0D0,0.0D0,Z1+ZL)
C      CALL LINE(IX1,IY1,IX2,IY2)
C      CALL CHARAC(IX2-IXX2,IY2+IXX1,CHARZ)
C      CALL LINE(IX2,IY2,IX2-IXX2,IY2-IXX4)
C      CALL LINE(IX2,IY2,IX2+IXX2,IY2-IXX4)
C      CALL LCOLOR(9)
C
      DO 1000 I=1,NT
         NT3=(I-1)*3
         DO 100 J=1,3
            II(J)=ITRI(NT3+J)
            XW(J)=VERTX(II(J))
            YW(J)=VERTY(II(J))
            IF (II(J).LT.0) THEN
               ZW(J)=0
            ELSE
               ZW(J)=U(II(J))
            END IF
  100    CONTINUE
         DO 200 J=1,3
            JJ=MOD(J,3)+1
C            CALL LINE(IGX(XW(J),YW(J),ZW(J)),IGY(XW(J),YW(J),ZW(J)),
C     +             IGX(XW(JJ),YW(JJ),ZW(JJ)),IGY(XW(JJ),YW(JJ),ZW(JJ)))
  200    CONTINUE
 1000 CONTINUE
      RETURN
      END
C
C Calculate the Value of the FEM Solution at a given Point.
      FUNCTION VALUE(X,Y,U,VERTX,VERTY,ITRI,NT,HX,HY)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER(NBPNO=-1300,NIPNO=500000,ITRINO=800000)
C      PARAMETER(NBPNO=-1000,NIPNO=30000,ITRINO=100000)
C      PARAMETER(NBPNO=-800,NIPNO=10000,ITRINO=60000)
      REAL*8 VERTX(NBPNO:NIPNO),VERTY(NBPNO:NIPNO),U(NIPNO),
     +       XW(3),YW(3),ZW(3)
      INTEGER*4 ITRI(ITRINO),II(3)
      Y0=Y-HY
      X0=X-HX
      X1=X+HX
      DO 200 I=1,NT
         NT3=(I-1)*3
         YW(1)=VERTY(ITRI(NT3+1))
         IF (YW(1).LT.Y0) GO TO 200
         XW(1)=VERTX(ITRI(NT3+1))
         IF (XW(1).LT.X0) GO TO 200
         IF (XW(1).GT.X1) GO TO 200
         DO 100 J=2,3
            XW(J)=VERTX(ITRI(NT3+J))
            YW(J)=VERTY(ITRI(NT3+J))
  100    CHI1=ORIENT(X,Y,XW(2),YW(2),XW(3),YW(3))
         IF (CHI1.LT.0.0D0) GO TO 200
         CHI2=ORIENT(X,Y,XW(3),YW(3),XW(1),YW(1))
         IF (CHI2.LT.0.0D0) GO TO 200
         CHI3=ORIENT(X,Y,XW(1),YW(1),XW(2),YW(2))
         IF (CHI3.LT.0.0D0) GO TO 200
         D=ORIENT(XW(1),YW(1),XW(2),YW(2),XW(3),YW(3))
         VALUE=(CHI1*U(ITRI(NT3+1))+CHI2*U(ITRI(NT3+2))
     +         +CHI3*U(ITRI(NT3+3)))/D
         WRITE(*,'(1H ,A17,I7)')'HITTED TRIANGLE :',I
         RETURN
  200 CONTINUE
      WRITE(*,*)'NO HITTED TRIANGLE'
      RETURN
      END
C
C Subroutine for 'PREPMV'
      FUNCTION ORIENT(X1,Y1,X2,Y2,X3,Y3)
      IMPLICIT REAL*8(A-H,O-Z)
      ORIENT=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
      RETURN
      END
C
C Function on the Right-Hand Side
      FUNCTION F(X,Y)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      F=5
      RETURN
      END