C*************************************************************
C     Test program for subroutine LSEMA by Stephen Kirkup            
C*************************************************************
C
C  Version 2. Revised test problems. Changes
C   in interface with LSEM3.
C  School Engineering, University of Central Lancashire 
C  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/LSEM3_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***************************************************************
C This program is a test for the subroutine LSEMA. The program computes
C  the solution to the Laplace equation 
C
C                  __ 2         
C                  \/    {\phi}   =  0  
C
C  exterior to a (set of) thin shell(s).
C
C Subroutine LSEMA accepts a description of the shell of the domain
C  and the position of the exterior points where the solution ({\phi}) 
C  is sought, the shell 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 two circular plates of
C  radius 1 (metre) and 0.1 apart. The solution to the problem with a
C  Dirichlet shell condition ({a}=1, {b}=0) with
C  {\phi}=f=1.0 on the upper plate and {\phi}=f=0.0 on the lower plate. 
C
C We would expect a simple gradient in potential between the two
C   plates and this is investigated in this test program.
C
C The plates are described by a set of NH=20 panels of equal length
C  along the generator. The shell solution points are the centres of the 
C  generator of the panels. 
C The solution is sought the points (0,0,0.025), (0,0,0.05), 
C  (0,0,0.075).

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

C The PARAMETER statement
C -----------------------
C There are four components in the PARAMETER statement.
C integer MAXNH   : The limit on the number of shell 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 LSEMA: Subroutine for solving the exterior Laplace
C  equation. (file LSEMA.FOR contains LSEMA 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 GLS2.FOR )
C file GEOM2D.FOR contains the set of relevant geometric subroutines 


C The program 

      PROGRAM LSEMAT
      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    MAXNH
      PARAMETER (MAXNH=32)
C   Limit on the number of vertices (equal to the number of panels)
      INTEGER    MAXNV
      PARAMETER (MAXNV=MAXNH)
C   Limit on the number of points exterior to the shell, where 
C    the solutions are sought
      INTEGER    MAXNPE
      PARAMETER (MAXNPE=6)

C  Geometrical description of the shell(s)
C   Number of panels and counter
      INTEGER    NH,IH
C   Number of central points (on H) and counter
      INTEGER    NHP,IHP
C   Number of vetices and counter
      INTEGER    NV,IV
C   Index of nodal coordinate for defining shell(s) standard unit is 
C    metres)
      REAL*8     VERTEX(MAXNV,2)
C   The two nodes that define each panel on the boundaries
      INTEGER    HELV(MAXNH,2)
C   The points exterior to the shell(s) 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 LSEMA.
C    HA(j) is assigned the value of {a} at the centre of the 
C     j-th panel.
      REAL*8     HA(MAXNH)
C    HB(j) is assigned the value of {b} at the centre of the 
C     j-th panel.
      REAL*8     HB(MAXNH)
C    HF(j) is assigned the value of f at the centre of the j-th panel.
      REAL*8     HF(MAXNH)
C    HA(j) is assigned the value of A at the centre of the 
C     j-th panel.
      REAL*8     HAA(MAXNH)
C    HB(j) is assigned the value of B at the centre of the 
C     j-th panel.
      REAL*8     HBB(MAXNH)
C    HF(j) is assigned the value of F at the centre of the j-th panel.
      REAL*8     HFF(MAXNH)

C The incident field
C  The incident potential on the shell(s)
      REAL*8     HIPHI(MAXNH)
C  The derivative of the incident potential on the shell(s)
      REAL*8     HIVEL(MAXNH)
C  The incident potential at the exterior points
      REAL*8     PEIPHI(MAXNH)

      
C  Validation and control parameters for LSEMA
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 shell.
      REAL*8     EGEOM

C Output from subroutine LSEMA
C  The average potential at the centres of the panels
      REAL*8     PHIAV(MAXNH)
C  The difference in potential at the centres of the panels
      REAL*8     PHIDIF(MAXNH)
C  The average derivative of the potential at the centres of the 
C   panels
      REAL*8     VELAV(MAXNH)
C  The difference in the derivative of the potential at the centres 
C   of the panels
      REAL*8     VELDIF(MAXNH)
      REAL*8     PEPHI(MAXNPE)

C  Working space 
      REAL*8     AMAT(2*MAXNH,2*MAXNH)
      REAL*8     BMAT(2*MAXNH,2*MAXNH)
      REAL*8     L_EH(MAXNPE,MAXNH)
      REAL*8     M_EH(MAXNPE,MAXNH)

C Further information from system solution GLS
      INTEGER*4  PERM(2*MAXNH)
      LOGICAL    XORY(2*MAXNH)
      REAL*8     C(2*MAXNH)

C Work Space
      REAL*8     WKSPC1(2*MAXNH)
      REAL*8     WKSPC2(2*MAXNH)
      REAL*8     WKSPC3(2*MAXNH)
C  Counter through the x,y coordinates
      INTEGER    ICOORD


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


C Describe the nodes and panels that make up the shell
C  :The plates are circles, each divided into 10 panels. VERTEX and 
C  : HELV are defined out from the centres of the plates so that the 
C  : normal to the plates is assumed to be upward (+z direction).
C  :Set up nodes
C  : Set NH, the number of panels
      NH=20
C  : Set NV, the number of vertices (equal to the number of panels)
      NV=22
C  : Set coordinates of the nodes
      DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,22)
     * / 0.000D0, 0.000D0,
     *   0.100D0, 0.000D0,
     *   0.200D0, 0.000D0,
     *   0.300D0, 0.000D0,
     *   0.400D0, 0.000D0,
     *   0.500D0, 0.000D0,
     *   0.600D0, 0.000D0,
     *   0.700D0, 0.000D0,
     *   0.800D0, 0.000D0,
     *   0.900D0, 0.000D0,
     *   1.000D0, 0.000D0,
     *   0.000D0, 0.100D0,
     *   0.100D0, 0.100D0,
     *   0.200D0, 0.100D0,
     *   0.300D0, 0.100D0,
     *   0.400D0, 0.100D0,
     *   0.500D0, 0.100D0,
     *   0.600D0, 0.100D0,
     *   0.700D0, 0.100D0,
     *   0.800D0, 0.100D0,
     *   0.900D0, 0.100D0,
     *   1.000D0, 0.100D0/

C  :Describe the panels that make up the two shells
C  : Set NH, the number of panels
      NH=20
C  : Set nodal indices that describe the panels of the shells.
C  :  The indices refer to the nodes in VERTEX. The order of the
C  :  nodes in HELV dictates that the normal is outward from the 
C  :  shell into the domain.
      DATA ((HELV(IH,ICOORD),ICOORD=1,2),IH=1,20)
     * /  1,  2,
     *    2,  3,
     *    3,  4,
     *    4,  5,
     *    5,  6,
     *    6,  7,
     *    7,  8,
     *    8,  9,
     *    9,  10,
     *   10,  11,
     *   12,  13,
     *   13,  14,
     *   14,  15,
     *   15,  16,
     *   16,  17,
     *   17,  18,
     *   18,  19,
     *   19,  20,
     *   20,  21,
     *   21,  22 /
       


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,     0.025D0,
     *     0.000D0,     0.050D0,
     *     0.000D0,     0.075D0/


C The number of points on the shell is equal to the number of 
C  panels
      NHP=NH
        
C  :In the test problem the plate 1 has                
C    {\phi}= 1
C  : and plate 2 has
C    {\phi}= 0
 
      DO 250 IHP=1,NHP/2
        HA(IHP)=1.0D0
        HB(IHP)=0.0D0
        HF(IHP)= 0.0D0
        HAA(IHP)=1.0D0
        HBB(IHP)=0.0D0
        HFF(IHP)= 0.0D0
250   CONTINUE     
      DO 260 IHP=NHP/2+1,NHP
        HA(IHP)=1.0D0
        HB(IHP)=0.0D0
        HF(IHP)= 0.0D0
        HAA(IHP)=1.0D0
        HBB(IHP)=0.0D0
        HFF(IHP)= 1.0D0
260   CONTINUE     

C There is no incident field
      DO 640 IHP=1,NHP
        HIPHI(IHP)=0.0D0
        HIVEL(IHP)=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 LSEMA
      LVALID=.TRUE.
C  :Set EGEOM
      EGEOM=1.0D-6


C   Set up particular a,b,f,A,B,F functions for this wavenumber
C    and type of shell condition

       
      CALL LSEMA(MAXNV,NV,VERTEX,MAXNH,NH,HELV,
     *                 MAXNPE,NPE,PEXT,
     *                 HA,HB,HF,HAA,HBB,HFF,
     *                 HIPHI,HIVEL,PEIPHI,
     *                 LSOL,LVALID,EGEOM,
     *                 PHIDIF,PHIAV,VELDIF,VELAV,PEPHI,
     *                 AMAT,BMAT,L_EH,M_EH,
     *                 PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)


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

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

      END

              