C***************************************************************
C     Test program for subroutine LEBEMA 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/LEBEMA_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 LEBEMA. The program computes
C  the solution to the Laplace equation 
C
C                  __ 2         
C                  \/    {\phi}   =  0  
C
C  exterior to a sphere by the boundary panel method where {\phi} is 
C  known as the potential.
C
C For the exterior problem, the domain lies exterior to a closed 
C  boundary S. The boundary condition may be Dirichlet, Robin or 
C  Neumann. It is assumed to have the following general form
C
C            {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q)
C    
C  where {\phi}(q) is the velocity potential at the point q on S, v(q) 
C  is the derivative of {\phi} with respect to the outward normal to S 
C  at q and {\alpha}, {\beta} and f are real-valued functions defined
C   on S. 
C
C Subroutine LEBEMA accepts a description of the boundary of the domain
C  and the position of the exterior points where the solution ({\phi}) 
C  is sought, the boundary condition and returns the solution ({\phi} 
C  and v) on S and the value of {\phi} at the exterior points.
C
C
C The test problems
C -----------------
C
C In this test the domain is the exterior to a sphere of side 1 
C  (metre). The solution to the problem with a Dirichlet boundary 
C  condition ({\alpha}=1, {\beta}=0) and with a Neumann boundary 
C  condition ({\alpha}=0, beta=1) are sought. 
C
C In the R-z plane, the boundary conditions are specified through 
C  taking the solution to be determined by
C
C              {\phi} = 1/R ,
C
C  which gives v=1 on the surface.
C
C The boundary is described by a set of NS=32 panels of equal size,
C  so that each side comprises eight panels. The boundary solution
C  points are the centres of the panels. 
C The *s show the exterior points at which the solution is sought;
C  the points (0,2), (1,1), (0,100).

C----------------------------------------------------------------------

C The PARAMETER statement
C -----------------------
C There are four components in the PARAMETER statement.
C integer MAXNS   : The limit on the number of boundary panels.
C integer MAXNV   : The limit on the number of vertices.
C integer MAXNPE  : The limit on the number of exterior points.


C External modules related to the package
C ---------------------------------------
C subroutine LEBEMA: Subroutine for solving the exterior Laplace
C  equation. (file LEBEMA.FOR contains LEBEMA and subordinate routines)
C subroutine L3ALC: Returns the individual discrete Laplace integral
C  operators. (file L3ALC.FOR contains L3ALC and subordinate routines)
C subroutine GLS: Solves a general linear system of equations.
C  (file GLS.FOR contains CGSL and subordinate routines)
C file GEOM2D.FOR contains the set of relevant geometric subroutines 


C The program 

      PROGRAM LEBEMAT
      IMPLICIT NONE

C VARIABLE DECLARATION
C --------------------

C  PARAMETERs for storing the limits on the dimension of arrays
C   Limit on the number of panels
      INTEGER    MAXNS
      PARAMETER (MAXNS=32)
C   Limit on the number of vertices (equal to the number of panels)
      INTEGER    MAXNV
      PARAMETER (MAXNV=MAXNS)
C   Limit on the number of points exterior to the boundary, where 
C    the solutions are sought
      INTEGER    MAXNPE
      PARAMETER (MAXNPE=6)

C  Constants
C   Real scalars: 0, 1, 2, pi
      REAL*8 ZERO,ONE,TWO,FOUR,PI

C  Geometrical description of the boundary(ies)
C   Number of panels and counter
      INTEGER    NS,IS
C   Number of central points (on S) and counter
      INTEGER    NSP,ISP
C   Number of vetices and counter
      INTEGER    NV,IV
C   Index of nodal coordinate for defining boundaries (standard unit is 
C    metres)
      REAL*8     VERTEX(MAXNV,2)
C   The two nodes that define each panel on the boundaries
      INTEGER    SELV(MAXNS,2)
C   The points exterior to the boundary(ies) where the solutions 
C    are sought and the directional vectors at those points.
C    [Only necessary if an exterior solution is sought.]
C    Number of exterior points and counter
      INTEGER    NPE,IPE
C    Coordinates of the exterior points
      REAL*8     PEXT(MAXNPE,2)


C   Data structures that are used to define each test problem in turn
C    and are input parameters to LEBEMA.
C    SALPHA(j) is assigned the value of {\alpha} at the centre of the 
C     j-th panel.
      REAL*8 SALPHA(MAXNS)
C    SBETA(j) is assigned the value of {\beta} at the centre of the 
C     j-th panel.
      REAL*8 SBETA(MAXNS)
C    SF(j) is assigned the value of f at the centre of the j-th panel.
      REAL*8 SF(MAXNS)

C  The incident field
C   The incident potential on the panels of the shell(s)
      REAL*8 SIPHI(MAXNS)
C   The derivative of the incident potential on the panels of the shell(s)
      REAL*8 SIVEL(MAXNS)
C   The incident potential at the domain points
      REAL*8 PEIPHI(MAXNS)

      
C  Validation and control parameters for LEBEMA
C   Switch for particular solution
      LOGICAL    LSOL
C   Validation switch
      LOGICAL    LVALID
C   The maximum absolute error in the parameters that describe the
C    geometry of the boundary.
      REAL*8     EGEOM

C Output from subroutine LEBEMA
C  The velocity potential (phi - the solution) at the centres of the 
C   panels
      REAL*8 SPHI(MAXNS)
C  The normal derivative of the velocity potential at the centres of the
C    panels
      REAL*8 SVEL(MAXNS)
C  The velocity potential (phi - the solution) at exterior points
      REAL*8 PEPHI(MAXNPE)

C Workspace for LEBEMA
      REAL*8 WKSPC1(MAXNS,MAXNS)
      REAL*8 WKSPC2(MAXNS,MAXNS)
      REAL*8 WKSPC3(MAXNPE,MAXNS)
      REAL*8 WKSPC4(MAXNPE,MAXNS)
      REAL*8 WKSPC5(MAXNS)
      REAL*8 WKSPC6(MAXNS)
      LOGICAL    WKSPC7(MAXNS)


C  Counter through the x,y coordinates
      INTEGER    ICOORD

C  The coordinates of the centres of the panels  
      REAL*8     SELCNT(MAXNS,2)

C  Other variables used in specifying the boundary condition
      REAL*8     Z,R

      REAL*8     EPS

C INITIALISATION
C --------------

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      FOUR=4.0D0
      PI=4.0D0*ATAN(ONE)

      EPS=1.0E-10


C Describe the nodes and panels that make up the boundary
C  :The circle that generates the sphere is divided into NS=18 uniform
C  : panels. VERTEX and SELV are defined anti-clockwise around the
C  : boundary so that the normal to the boundary is assumed to be 
C  : outward
C  :Set up nodes
C  : Set NS, the number of panels
      NS=18
C  : Set NV, the number of vertices (equal to the number of panels)
      NV=NS+1
C  : Set coordinates of the nodes
      DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,19)
     * / 0.000D0, 1.000D0,
     *   0.174D0, 0.985D0,
     *   0.342D0, 0.940D0,
     *   0.500D0, 0.866D0,
     *   0.643D0, 0.766D0,
     *   0.766D0, 0.643D0,
     *   0.866D0, 0.500D0,
     *   0.940D0, 0.342D0,
     *   0.985D0, 0.174D0,
     *   1.000D0, 0.000D0,
     *   0.985D0,-0.174D0,
     *   0.940D0,-0.342D0,
     *   0.866D0,-0.500D0,
     *   0.766D0,-0.643D0,
     *   0.643D0,-0.766D0,
     *   0.500D0,-0.866D0,
     *   0.342D0,-0.940D0,
     *   0.174D0,-0.985D0,
     *   0.000D0,-1.000D0 /

C  :Describe the panels that make up the two boundarys
C  : Set NS, the number of panels
      NS=18
C  : Set nodal indices that describe the panels of the boundarys.
C  :  The indices refer to the nodes in VERTEX. The order of the
C  :  nodes in SELV dictates that the normal is outward from the 
C  :  boundary into the domain.
      DATA ((SELV(IS,ICOORD),ICOORD=1,2),IS=1,18)
     * /  1,  2,
     *    2,  3,
     *    3,  4,
     *    4,  5,
     *    5,  6,
     *    6,  7,
     *    7,  8,
     *    8,  9,
     *    9,  10,
     *   10,  11,
     *   11,  12,
     *   12,  13,
     *   13,  14,
     *   14,  15,
     *   15,  16,
     *   16,  17,
     *   17,  18,
     *   18,  19 /
       
C Set the centres of the panels, the central points
      DO IS=1,NS
        SELCNT(IS,1)=(VERTEX(SELV(IS,1),1)
     *   +VERTEX(SELV(IS,2),1))/TWO
        SELCNT(IS,2)=(VERTEX(SELV(IS,1),2)
     *   +VERTEX(SELV(IS,2),2))/TWO
      END DO


C Set the points in the domain where the solutions are sought, PEXT. 
      NPE=3
      DATA ((PEXT(IPE,ICOORD),ICOORD=1,2),IPE=1,3)
     *  /  0.000D0,     2.000D0,
     *     1.000D0,     1.000D0,
     *     0.000D0,    100.000D0/


C The number of points on the boundary is equal to the number of 
C  panels
      NSP=NS
        


C  :The test problem is an pulsating sphere so that                
C    {\phi}= 1
C  :This gives the following v on S 
C     v = -1
 
      DO 250 ISP=1,NSP/2
        R=SELCNT(ISP,1)
        Z=SELCNT(ISP,2)
        SALPHA(ISP)=1.0D0
        SBETA(ISP)=0.0D0
        SF(ISP)= 1
250   CONTINUE     
      DO 260 ISP=NSP/2+1,NSP
        R=SELCNT(ISP,1)
        Z=SELCNT(ISP,2)
        SALPHA(ISP)=0.0D0
        SBETA(ISP)=1.0D0
        SF(ISP)= -1
260   CONTINUE     
      DO 640 ISP=1,NSP
        SIPHI(ISP)=0.0D0
        SIVEL(ISP)=0.0D0
640   CONTINUE
      DO 650 IPE=1,NPE
        PEIPHI(IPE)=0.0D0
650   CONTINUE


C Set up validation and control parameters
C  :Switch for particular solution
      LSOL=.TRUE.
C  :Switch on the validation of LEBEMA
      LVALID=.TRUE.
C  :Set EGEOM
      EGEOM=1.0D-6


C   Set up particular alpha and beta functions for this wavenumber
C    and type of boundary condition

       
          CALL LEBEMA(MAXNV,NV,VERTEX,MAXNS,NS,SELV,
     *               MAXNPE,NPE,PEXT,
     *               SALPHA,SBETA,SF,SIPHI,SIVEL,PEIPHI,
     *               LSOL,LVALID,EGEOM,
     *               SPHI,SVEL,PEPHI,
     *               WKSPC1,WKSPC2,WKSPC3,WKSPC4,
     *               WKSPC5,WKSPC6,WKSPC7)


C Output the solutions
C  Open file for the output data
      OPEN(UNIT=20,FILE='LEBEMA.OUT')

      WRITE(20,*) 'LEBEMA: Computed solutions in the domain'
      WRITE(20,*)
      WRITE(20,*) '   P(1)      P(2)      Computed    Exact'
      DO 800 IPE=1,NPE
        WRITE(20,999) PEXT(IPE,1),PEXT(IPE,2),PEPHI(IPE),
     *   1.0/SQRT(PEXT(IPE,1)**2+PEXT(IPE,2)**2)
800   CONTINUE
999   FORMAT(4F10.4)
      CLOSE(20)

      END

                                                           