C***************************************************************
C    Test program for subroutine L2LC by Stephen Kirkup               
C***************************************************************
C 
C  Copyright 2001- Stephen Kirkup
C  School of Computing Engineering and Physical Sciences
C  University of Central Lancashire - www.uclan.ac.uk 
C  smkirkup@uclan.ac.uk
C  http://www.researchgate.net/profile/Stephen_Kirkup
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/L2LC_T.FOR 
C
C  Issued under the GNU General Public License 2007, see gpl.txt
C
C  Part of the the author's open source BEM packages. 
C  All codes and manuals can be downloaded from 
C  www.boundary-element-method.com
C
C***************************************************************

C This program is a test for the subroutine L2LC. The program computes
C  the solution to the Laplace problem interior In order to use L2LC, 
C  the circle is approximated by a regular polygon with each side being
C  one element. The boundary functions are approximated by a constant 
C  at the centre of each element.
C Simple changes in the program allow the set up of the following data
C  (1) the choice of quadrature rule,
C  (2) the radius of the circle,
C  (3) the number of elements (sides on the approximating polygon),
C  (4) the position of one source point and one sink point which 
C       determine the boundary condition and the exact solution,
C  (5) the points in the interior where the solution is sought.
C During execution the program gives the solution at the collocation
C  points (the points at the centre of each element) and the solution
C  at the selected interior points. The program also give the exact
C  solution at the same points so that computed and exact solutions
C  may be compared.

C The PARAMETER statement
C =======================
C There are four components in the PARAMETER statement.
C integer LIMN   : The limit on the number of elements.
C integer LIMNQ  : The limit on the number of quadrature points allowed
C                   in the quadrature rules.
C integer LIMNPI : The limit on the number of interior points.

C External modules related to the test program
C ============================================
C Subroutine GLRULE: This generates Gauss-Legendre quadrature rules.
C real function FNLOG: Returns the natural logarithm.
C real function FNSQRT: Returns the square root.

C External modules related to the package
C =======================================
C Subroutine L2LC: This computes the discrete Laplace integral
C  operators.

      PROGRAM     IBEM2
C Declaration of variables
C  Limits on the size of the arrays
C   Limit on the number of elements
      INTEGER     LIMN
C   Limit on the number of quadrature points in any quadrature rule
      INTEGER     LIMNQ
C   Limit on the number of interior points where a solution is sought
      INTEGER     LIMNPI
C  Set PARAMETERs
      PARAMETER (LIMN=64,LIMNQ=16,LIMNPI=4)
C  Constants
      REAL*8     ZERO,HALF,ONE,TWO,FOUR,PI,TWOPI
      REAL*8 CIMAG
C  Number of elements
      INTEGER    N
C  Geometrical information
      INTEGER    NPEXT
      REAL*8     P(2),NORMP(2),Q(2),QA(2),QB(2),VECP(2)
      REAL*8     NORMQ(2),PINT(LIMNPI,2)
      REAL*8     SEGANG,ANGL,ANGU,RAD,ANGP
      REAL*8     ANG
      LOGICAL    LPONEL
C  Error flag used for checking the set-up information
      LOGICAL    LERROR
C  Quadrature rules
C   General quadrature rule
C    Number of points
      INTEGER    NQ
C    Points and weights
      REAL*8     AQ(LIMNQ),WQ(LIMNQ)
C   Particular quadrature rule 
C    Number of points
      INTEGER    NQT
C    Points and weights
      REAL*8     AQT(LIMNQ),WQT(LIMNQ)
C  Validation and control variables in L2LC
      LOGICAL    LVALID,LFAIL
      REAL*8     EGEOM,EQRULE
      LOGICAL    LL,LM,LMT,LN
C  Storage of discrete integral operators in L2LC
      REAL*8     DISL,DISM,DISMT,DISN
C  Matrices corresponding to the Laplace integral operators
      REAL*8     L(LIMN,LIMN),M(LIMN,LIMN)
C  Composite matrices 
      REAL*8     A(LIMN,LIMN),B(LIMN,LIMN)
C  Exact Dirichlet boundary condition, exact solution and computed 
C   solution
      REAL*8     PHI(LIMN),VEL(LIMN),APPVEL(LIMN)
C  Solution at the interior point
      REAL*8     PHIPT
C  Temporary variable used in computing the solutions
      REAL*8     APHI(LIMN),SUM
C  Counters and indices
      INTEGER    I,J,IQ,IPI
      
C Set up constants
      ZERO=0.0D0
      HALF=0.5D0
      ONE=1.0D0
      TWO=2.0D0
      FOUR=4.0D0
      CIMAG=CMPLX(ZERO,ONE)
      PI=FOUR*ATAN(ONE)
      TWOPI=TWO*PI

C Set up control and tolerances and open file for storing the possible
C  error messages
      LVALID=.TRUE.
      EGEOM=1D-08
      EQRULE=1D-08
      OPEN(UNIT=10,FILE='L2LC.ERR',STATUS='UNKNOWN')

C Set up standard Gaussian Quadrature rule
      NQ=8
      CALL GLRULE(NQ,ZERO,ONE,WQ,AQ)

C Set radius of circle (RAD)
      RAD=1.0D0

C Set the number of elements (N)
      N=64

C Set interior points, where an approximation the the solution phi is 
C  sought
      NPEXT=1
      PINT(1,1)=0.0D0
      PINT(1,2)=0.5D0

C Check that the set-up data is satisfactory
      LERROR=.FALSE.
      IF (NQ.GT.LIMNQ) THEN
        WRITE(*,*) 'ERROR - NQ must be less than or equal to LIMNQ'
        LERROR=.TRUE.
      END IF
      IF (RAD.LE.ZERO) THEN
        WRITE(*,*) 'ERROR - RAD must be greater than 0'
        LERROR=.TRUE.
      END IF
      IF (N.GT.LIMN) THEN
        WRITE(*,*) 'ERROR - N must be less than or equal to LIMN'
        LERROR=.TRUE.
      END IF
      IF (NPEXT.GT.LIMNPI) THEN
        WRITE(*,*) 'ERROR - NPEXT must be less than or equal to LIMNPI'
        LERROR=.TRUE.
      END IF
      DO 50 IPI=1,NPEXT
        IF (PINT(IPI,1)*PINT(IPI,1)+PINT(IPI,2)*PINT(IPI,2).GT.RAD*RAD)
     *    THEN
          WRITE(*,*) 'ERROR - PINT must be interior to the circle'
          LERROR=.TRUE.
        END IF
50    CONTINUE
      IF (LERROR) STOP
      

C Output description of test the problem
      WRITE(*,*) 'TEST FOR SUBROUTINE L2LC'
      WRITE(*,*) '========================='
      WRITE(*,*) 'Test problem is that of a circle centred at the'
      WRITE(*,*) 'origin with a Dirichlet boundary condition'
      WRITE(*,*) 'The radius of the circle is ',RAD
      WRITE(*,*) 'Number of boundary elements is',N
      

C  The point P is defined by setting the angle to -SEGANG/2 then adding
C   SEGANG systematically to give the angle each point makes with the 
C   y-axis at the origin.
        SEGANG=TWOPI/DBLE(N)
        ANGP=-SEGANG/TWO
C  Loop(I) through each collocation point 
        DO 100 I=1,N

C   Set P(1..2) (the collocation point) and NORMP(1..2) the normal at 
C    P(1..2).
          ANGP=ANGP+SEGANG
          P(1)=RAD*COS(SEGANG/TWO)*SIN(ANGP)
          P(2)=RAD*COS(SEGANG/TWO)*COS(ANGP)
          NORMP(1)=SIN(ANGP)
          NORMP(2)=COS(ANGP)

          ANGU=0.0D0
C   Loop(J) through the elements
          DO 110 J=1,N

C    Set QA(1..2) and QB(1..2) the vertices of the element. Since the
C     integral that is applied assumes the normal to be  outward,
C     the vertices rotate clockwise around the boundary. This is
C     so that the normal is implicitly defined in L2LC as being 
C     outward from the boundary of the circle. 
            ANGL=ANGU
            ANGU=ANGL+SEGANG
            QA(1)=RAD*SIN(ANGL)
            QA(2)=RAD*COS(ANGL)
            QB(1)=RAD*SIN(ANGU)
            QB(2)=RAD*COS(ANGU)

C    Set LPONEL
            IF (I.EQ.J) THEN
              LPONEL=.TRUE.
            ELSE
              LPONEL=.FALSE.
            END IF

C    Select the quadrature rule and store in NQT,AQT,WQT. If LPONEL is
C     false then use the existing quadrature rule NQ,AQ,WQ. If LPONEL
C     is true then split the NQ,AQ,WQ quadrature rule into two at the 
C     centre. This deals with the dicontinuity at the centre.
            IF (LPONEL) THEN
              NQT=2*NQ
              DO 120 IQ=1,NQ
                WQT(IQ)=WQ(IQ)/TWO
                WQT(NQ+IQ)=WQ(IQ)/TWO
                AQT(IQ)=AQ(IQ)/TWO
                AQT(NQ+IQ)=(ONE+AQ(IQ))/TWO
120           CONTINUE
              ELSE
              NQT=NQ
              DO 130 IQ=1,NQ
                WQT(IQ)=WQ(IQ)
                AQT(IQ)=AQ(IQ)
130           CONTINUE
            ENDIF

C    Only integral operators L and M are required
            LL=.TRUE.
            LM=.TRUE.
            LMT=.FALSE.
            LN=.FALSE.

C    Call of L2LC routine to compute [Lk], [Mk]
          CALL L2LC(P,NORMP,QA,QB,LPONEL,
     *     LIMNQ,NQT,AQT,WQT,
     *     LVALID,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)

C    Storage of results in matrices
          L(I,J)=DISL
          M(I,J)=DISM

C   End loop(J) through elements
110     CONTINUE

C  End loop(I) through collocation points
100   CONTINUE

C  Set up test Dirichlet boundary condition (PHI) and the exact solution
C   on the boundary (VEL). Here \phi = y which has the exact solution
C   v=y
      ANG=-SEGANG/TWO
      DO 140 I=1,N
        ANG=ANG+SEGANG
        Q(1)=RAD*SIN(ANG)
        Q(2)=RAD*COS(ANG)
        NORMQ(1)=SIN(ANG)
        NORMQ(2)=COS(ANG)
        PHI(I)=Q(2)
        VEL(I)=Q(2)
140   CONTINUE

C  Construct composite matrices A and B for the linear system
      DO 150 I=1,N
        DO 160 J=1,N
          A(I,J)=M(I,J)
          B(I,J)=L(I,J)
160     CONTINUE
        A(I,I)=A(I,I)+HALF
150   CONTINUE

C  Solve linear system B APPVEL = A PHI to give APPVEL, the 
C   approximation to PHI.
      DO 170 I=1,N
        APHI(I)=ZERO
        DO 180 J=1,N
          APHI(I)=APHI(I)+A(I,J)*PHI(J)
180     CONTINUE
170   CONTINUE

C   Call the subroutine for solving the linear system of equations
      CALL LINSL(LIMN,N,B,APPVEL,APHI)

C  Output of solution on boundary
999   FORMAT(2F12.6,'   ',2F12.6)
      WRITE(*,*) '      EXACT PHI ON S            APPROX PHI ON S'
      DO 190 I=1,N
        WRITE(*,999) VEL(I), APPVEL(I)
190   CONTINUE

C  Loop(IPI) through the points in the interior
      DO 400 IPI=1,NPEXT
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)

C   Output the coordinates of the point in the interior
        WRITE(*,*) 'Point in the interior is ',P

C   Compute the exact solution PHIPT at the interior point
C    Exact solution is \phi = y
        PHIPT=P(2)

C   Initialise SUM, this ultimately stores the value of the summation
C    that approximates to phi at the interior point
        SUM=ZERO
        ANGU=0.0D0

C   Loop(J) through the elements 
        DO 200 J=1,N

C    Set QA(1..2) and QB(1..2) the vertices of the element.
          ANGL=ANGU
          ANGU=ANGL+SEGANG
          QA(1)=RAD*SIN(ANGL)
          QA(2)=RAD*COS(ANGL)
          QB(1)=RAD*SIN(ANGU)
          QB(2)=RAD*COS(ANGU)

C    Set up the quadrature rule
          NQT=NQ
          DO 210 IQ=1,NQ
            WQT(IQ)=WQ(IQ)
            AQT(IQ)=AQ(IQ)
210       CONTINUE

C    For secondary stage of method Lk and Mk are required.
          LL=.TRUE.
          LM=.TRUE.
          LMT=.FALSE.
          LN=.FALSE.

C    Since P is in the interior then LPONEL is always false
          LPONEL=.FALSE.

C    Value of VECP is not relevant

C    Call of L2LC routine to compute [Lk], [Mk]
          CALL L2LC(P,VECP,QA,QB,LPONEL,
     *     LIMNQ,NQT,AQT,WQT,
     *     LVALID,EGEOM,EQRULE,LFAIL,
     *     LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)

C    Accumulate SUM
          SUM=SUM-DISM*PHI(J)+DISL*APPVEL(J)

C   End loop(J) through elements
200     CONTINUE

C   Ouput the solution at the interior point
        WRITE(*,*) 'EXACT  PHI AT POINT = ',PHIPT
        WRITE(*,*) 'APPROX PHI AT POINT = ',SUM

C  End loop(IPI) through interior points
400   CONTINUE


      CLOSE(10)

      END


C Subroutine LINSL
C ================

C This subroutine implements a Gaussian elimination procedure with 
C partial pivotting in order to solve the linear system of equations
C A X = B. AMAT(MAXN,MAXN) stores the NxN maxtrix on input but is altered
C on output. BVEC(MAXN) stores the N vector on input and XVEC(MAXN) the 
C solution on output.



        SUBROUTINE LINSL(MAXN,N,AMAT,XVEC,BVEC)
        IMPLICIT REAL*8 (A-H,O-Z)
        INTEGER    MAXN,N
        REAL*8 AMAT(MAXN,MAXN),BVEC(MAXN),XVEC(MAXN)
        REAL*8 TEMP,CSUM,RATIO
        REAL*8 EPS,AA,COLMAX
        INTEGER II,IIMAX,J,I
        EPS=1.0D-20

        DO 10 I=1,N-1
          COLMAX=0
          DO 20 II=I,N
            AA=ABS(AMAT(II,I))
            IF (AA.GT.COLMAX) THEN
              COLMAX=AA
              IIMAX=II
            ENDIF
20        CONTINUE
          DO 30 J=I,N
            TEMP=AMAT(I,J)
            AMAT(I,J)=AMAT(IIMAX,J)
            AMAT(IIMAX,J)=TEMP
30        CONTINUE
          DO 40 II=I+1,N
            IF (ABS(AMAT(I,I)).LT.EPS) THEN
              WRITE(*,*) 'ERROR(LINSL) - Input matrix A is singular'
              GOTO 998
            END IF
            RATIO=AMAT(II,I)/AMAT(I,I)
            DO 50 J=I+1,N
              AMAT(II,J)=AMAT(II,J)-RATIO*AMAT(I,J)
50          CONTINUE
            BVEC(II)=BVEC(II)-RATIO*BVEC(I)
40        CONTINUE
10      CONTINUE
        XVEC(N)=BVEC(N)/AMAT(N,N)
        DO 60 I=N-1,1,-1
          CSUM=BVEC(I)
          DO 70 J=I+1,N
            CSUM=CSUM-AMAT(I,J)*XVEC(J)
70        CONTINUE
          XVEC(I)=CSUM/AMAT(I,I)
60      CONTINUE

998     CONTINUE
        END





C Subroutine GLRULE: This generates Gauss-Legendre quadrature rules.
C This subroutine supplies the weights and abscissae for an N point
C Gaussian quadrature rule on the interval [A,B].

C Input parameters
C integer N   : The number of Gaussian points. N must satisfy the condions
C  required of its corresponding value in the NAG routine D01BBF.
C real*8  A   : The lower limit of the quadrature rule.
C real*8  B   : The upper limit of the quadrature rule.
C Output parameters
C real*8  WTS(N) : The Gaussian quadrature weights.
C real*8  PTS(N) : The Gaussian quadraure points.

C External modules
C D01BBF : From the NAG library.
C D01BAZ : From the NAG library.


      SUBROUTINE GLRULE(N,A,B,WTS,PTS)
      INTEGER N
      REAL*8 A,B,WTS(N),PTS(N)
      INTEGER IFAIL,ITYPE
      EXTERNAL D01BAZ
      ITYPE=1
      IFAIL=0
      CALL D01BBF(D01BAZ,A,B,ITYPE,N,WTS,PTS,IFAIL)
      END


C real function FNLOG: Returns the natural logarithm of X.
      REAL*8 FUNCTION FNLOG(X)
      REAL*8 X
      FNLOG=LOG(X)
      END


C real function FNSQRT: Returns the square root of X.
      REAL*8 FUNCTION FNSQRT(X)
      REAL*8 X
      FNSQRT=SQRT(X)
      END

                                                                             