C **************************************************************
C *  Test program for subroutine L3ALC 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/L3ALC.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 L3ALC. The program computes
C  the solution of the Laplace problem exterior to a sphere centred
C  at the origin and with an axisymmetric solution via the integral
C  equation (M+I/2)\phi = L v. The boundary condition and solution are 
C  assumed to be independent of theta (ie axisymmetric).
C In order to use L3ALC, the sphere is approximated by a set of conical
C  elements. Each element is decribed by a straight line on the generator 
C  (R-z plane) swept through 2 pi (in the theta direction). The boundary 
C  functions are approximated by a constant 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 sphere,
C  (3) the number of elements,
C  (4) the set of wavenumbers,
C  (5) the position of one source point and one sink point which 
C       determine the boundary condition and the exact solution,
C  (6) the points in the exterior 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 at
C  the selected exterior points. The program also give the exact
C  solution at the same points so that computed and exact solutions may
C  be compared.

C The PARAMETER statement
C =======================
C There are five components in the PARAMETER statements.
C integer LIMN  : The limit on the number of elements.
C integer LIMNGQ: The limit on the number of quadrature points allowed
C                  in the quadrature rules in the generator direction.
C integer LIMNTQ: The limit on the number of quadrature points allowed
C                  in the quadrature rules in the theta direction.
C integer LIMNPI: The limit on the number of interior points.

C External modules related to the test program
C ============================================
C Subroutine LINSL: This computes the solution to a linear system of
C  equations
C Subroutine GLRULE: This generates the Gauss-Legendre quadrature rules.
C real function FNSQRT: Returns the square root.

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

       PROGRAM    IBEM3A
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 generator-direction quadrature points in any
C    quadrature rule
       INTEGER    LIMNGQ
C   Limit on the number of theta-direction quadrature points in any
C    quadrature rule
       INTEGER    LIMNTQ
C   Limit on the number of interior points where a solution is sought
       INTEGER    LIMNPi
C  Set PARAMETERs 
       PARAMETER (LIMN=16,LIMNGQ=32,LIMNPI=2)
       PARAMETER (LIMNTQ=LIMN*LIMNGQ)
C  Constants
       REAL*8     ZERO,HALF,ONE,TWO,FOUR,EIGHT,PI,ROOT2
       REAL*8 CIMAG
C  Number of elements
       INTEGER    N
C  Geometrical information
       INTEGER    NPINT
       REAL*8     P(2),NORMP(2),Q(2),QA(2),QB(2)
       REAL*8     PINT(LIMNPI,2)
       REAL*8     SEGANG,ANGL,ANGU,ANGM,RADMID,GLEN,SGLEN,CIRMID,TDIV
       REAL*8     RAD,ANGP,VECP(2)
       LOGICAL    LPONEL
C  Error flag 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(LIMNGQ),WQ(LIMNGQ)
C   Generator-direction quadrature rule
C    Number of points
       INTEGER    NGQ
C    Points and weights
       REAL*8     AGQ(LIMNGQ),WGQ(LIMNGQ)
C    Number subdivisions of composite rule
       INTEGER    NDIV
C   Theta-direction quadrature rule
C    Number of points
       INTEGER    NTQ
C    Points and weights
       REAL*8     ATQ(LIMNTQ),WTQ(LIMNTQ)
C  Validation and control variables in L3ALC
       LOGICAL    LVALID,LFAIL
       REAL*8     EGEOM,EQRULE
       LOGICAL    LL,LM,LMT,LN
C  Storage of discrete integral operators in L3ALC
       REAL*8 DISL,DISM,DISMT,DISN
C  Matrices corresponding to the Laplace integral operators
       REAL*8 M(LIMN,LIMN),L(LIMN,LIMN)
C  Composite matrices used in the Burton and Miller equation
       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 exterior point
       REAL*8 PHIPT
C  Temporary variables used in computing the solutions
       REAL*8 APHI(LIMN),SUM
C  Counters and indices
       INTEGER*4  I,IDIV,IQ,J,IPE
C  Work space parameter for L3ALC
       REAL*8     WKSPCE(2*LIMNTQ+LIMNGQ)

C Set up constants
       ZERO=0.0D0
       HALF=0.5D0
       ONE=1.0D0
       TWO=2.0D0
       FOUR=4.0D0
       EIGHT=8.0D0
       CIMAG=CMPLX(ZERO,ONE)
       PI=FOUR*ATAN(ONE)
       ROOT2=SQRT(TWO)

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

C Set up standard Gaussian Quadrature rule. This quadrature rule is 
C  used directly along the generator for most integrals. In the
C  case when the point P lies at the centre of the element a composite
C  rule consisting of two sets of the standard rule is used on the
C  generator. For the theta direction rule a composite rule is formed
C  from the standard rule so that the density of quadrature points in
C  the theta direction is approximately the same as the density of
C  points in the generator direction.
       NQ=8
       CALL GLRULE(NQ,ZERO,ONE,WQ,AQ)

C Set radius of the sphere (RAD)
       RAD=1.0D0

C Set number of elements (N)
       N=16
        

C Set exterior points (R,z coordinates), where an approximation to the
C  solution phi is sought
       NPINT=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.LIMNGQ) THEN
         WRITE(*,*) 'ERROR - NQ must be less than or equal to LIMNGQ'
         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 (NPINT.GT.LIMNPI) THEN
         WRITE(*,*) 'ERROR - NPINT must be less than or equal to LIMNPI'
         LERROR=.TRUE.
       END IF
       DO 50 IPE=1,NPINT
         IF (PINT(IPE,1)*PINT(IPE,1)+PINT(IPE,2)*PINT(IPE,2).GT.RAD*RAD)
     *     THEN
           WRITE(*,*) 'ERROR - PINT must be interior to the circle'
           LERROR=.TRUE.
         END IF
50     CONTINUE
       IF (NQ*N.GT.LIMNTQ) THEN
         WRITE(*,*) 'ERROR - LIMNTQ should be at least NQ*N'
         LERROR=.TRUE.
       END IF
       IF (LERROR) STOP

C Output description of test the problem
       WRITE(*,*) 'TEST FOR SUBROUTINE L3ALC'
       WRITE(*,*) '========================='
       WRITE(*,*) 'Test problem is that of a sphere centred at the'
       WRITE(*,*) 'origin with a Dirichlet boundary condition'
       WRITE(*,*) 'The radius of the sphere 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 that each point makes with
C   the z-axis at the origin.
         SEGANG=PI/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 L3ALC as being 
C     outward from the boundary of the sphere. 
             ANGL=ANGU
             ANGU=ANGL+SEGANG
             ANGM=(ANGU+ANGL)/TWO
             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 generator-direction quadrature rule and store in NQT,
C     AQT,WQT. If LPONEL is false then use the existing quadrature rule 
C     NQ,AQ,WQ. If LPONEL is true then split the NQ,AQ,WQ quadrature 
C     rule into two at the centre. This deals with the dicontinuity at 
C     the centre.
             IF (LPONEL) THEN
               NGQ=2*NQ
               DO 120 IQ=1,NQ
                 WGQ(NQ+1-IQ)=WQ(IQ)/TWO
                 WGQ(2*NQ+1-IQ)=WQ(IQ)/TWO
                 AGQ(NQ+1-IQ)=AQ(IQ)/TWO
                 AGQ(2*NQ+1-IQ)=(ONE+AQ(IQ))/TWO
120            CONTINUE
             ELSE
               NGQ=NQ
               DO 130 IQ=1,NQ
                 WGQ(NQ+1-IQ)=WQ(IQ)
                 AGQ(NQ+1-IQ)=AQ(IQ)
130            CONTINUE
             ENDIF

C Quadrature rule in the theta direction is constructed out of individual
C Gauss rules so that the length of each is approximately equal to the
C length of the element at the generator.
             RADMID=(QA(1)+QB(1))/TWO
             SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *        (QA(2)-QB(2))*(QA(2)-QB(2))
             GLEN=SQRT(SGLEN)
             CIRMID=PI*RADMID
             NDIV=1+CIRMID/GLEN
             TDIV=ONE/DBLE(NDIV)
             NTQ=NDIV*NQ
             DO 140 IDIV=1,NDIV
               DO 150 IQ=1,NQ
                 WTQ((IDIV-1)*NQ+IQ)=WQ(IQ)/DBLE(NDIV)
                 ATQ((IDIV-1)*NQ+IQ)=AQ(IQ)/DBLE(NDIV)+TDIV*DBLE(IDIV-1)
150            CONTINUE
140          CONTINUE

C    For primary stage of method all integral operators are required
C     in general (for the Burton and Miller equation)
             LL=.TRUE.
             LM=.TRUE.
             LMT=.FALSE.
             LN=.FALSE.

C    Call of L3ALC routine to compute [Lk], [Mk]
             CALL L3ALC(P,NORMP,QA,QB,LPONEL,
     *        LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ,
     *        LVALID,EGEOM,EQRULE,LFAIL,
     *        LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *        WKSPCE)

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

C   End loop(J) through elements
110        CONTINUE

C  End loop(I) through collocation points
100      CONTINUE

C  Set up the Neumann boundary condition and the exact solution
C  Here \phi = z which has the solution v = z.
         ANGU=0.0D0
         DO 160 J=1,N
           ANGL=ANGU
           ANGU=ANGL+SEGANG
           ANGM=(ANGU+ANGL)/TWO
           Q(1)=RAD*SIN(ANGM)
           Q(2)=RAD*COS(ANGM)
           PHI(J)=Q(2)
           VEL(J)=Q(2)
160      CONTINUE

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

C  Solve linear system A PHI = B APPVEL to give APPVEL, the 
C   approximation to VEL.
         DO 190 I=1,N
           APHI(I)=ZERO
           DO 200 J=1,N
             APHI(I)=APHI(I)+A(I,J)*PHI(J)
200        CONTINUE
190      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, APPROX  VEL ON S'
         DO 210 I=1,N
           WRITE(*,999) VEL(I),APPVEL(I)
210      CONTINUE


C  Loop(IPE) through the points in the interior
         DO 300 IPE=1,NPINT
           P(1)=PINT(IPE,1)
           P(2)=PINT(IPE,2)

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

C   Compute the exact solution PHIPT at the exterior point
C   Here \phi = z
           PHIPT=P(2)

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

C   Loop(J) through the elements 
           DO 320 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
C     Generator direction
             NGQ=NQ
             DO 330 IQ=1,NQ
               WGQ(IQ)=WQ(NQ+1-IQ)
               AGQ(IQ)=AQ(NQ+1-IQ)
330          CONTINUE

C     Quadrature rule in the theta direction is constructed out of
C      individual Gauss rules so that the length of each is 
C      approximately equal to the length of the element at the 
C      generator.
             RADMID=(QA(1)+QB(1))/TWO
             SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+
     *        (QA(2)-QB(2))*(QA(2)-QB(2))
             GLEN=SQRT(SGLEN)
             CIRMID=PI*RADMID
             NDIV=1+CIRMID/GLEN
             TDIV=ONE/DBLE(NDIV)
             NTQ=NDIV*NQ
             DO 340 IDIV=1,NDIV
               DO 350 IQ=1,NQ
                 WTQ((IDIV-1)*NQ+IQ)=WQ(IQ)/DBLE(NDIV)
                 ATQ((IDIV-1)*NQ+IQ)=AQ(IQ)/DBLE(NDIV)+TDIV*DBLE(IDIV-1)
350            CONTINUE
340          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 exterior then LPONEL is always false
             LPONEL=.FALSE.

C    Value of VECP is not relevant

C    Call of L3ALC routine to compute [Lk], [Mk]
             CALL L3ALC(P,VECP,QA,QB,LPONEL,
     *        LIMNGQ,NGQ,AGQ,WGQ,LIMNTQ,NTQ,ATQ,WTQ,
     *        LVALID,EGEOM,EQRULE,LFAIL,
     *        LL,LM,LMT,LN,DISL,DISM,DISMT,DISN,
     *        WKSPCE)

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

C   End loop(J) through elements
320        CONTINUE

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

C  End loop(IPE) through the exterior points
300      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 the 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 routines
C D01BBF : From the NAG library.
C D01BAZ : From the NAG library.


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


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



                                       