C***************************************************************
C     Test program for subroutine LIBEM2 by Stephen Kirkup           
C***************************************************************
C 
C  Version 2 (Sept 2015) (tidy up of previous version Jul 2015)
C  (previous version 2001) Stephen Kirkup
C  School of Engineering 
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/LIBEM2_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 LIBEM2. The program computes
C  the solution to an Laplace problem interior to a square 
C  by the boundary panel method.
C
C Background
C ----------
C
C  The test problem solves the Laplace equation
C
C                  __ 2                
C                  \/    {\phi}   =  0  
C
C  where {\phi} is known as the potential.
C
C For the interior problem, the domain lies interior 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 complex-values functions defined
C   on S. 
C
C Subroutine LIBEM2 accepts a description of the 
C  boundary of the domain and the position of the interior points
C  where the solution ({\phi}) is sought, the boundary condition and
C  returns the solution ({\phi} and v) on S and the value of {\phi}
C  at the interior points.
C

C The test problems
C -----------------
C
C In this test the domain is a square of side 0.1 (metre) with a 
C  Dirichlet boundary condition ({\alpha}=1, {\beta}=0).
C Assuming the square has vertices (0,0), (0,0.1), (0.1,0.1) and (0.1,0)
C  in the x-y plane, the boundary conditions are specified through 
C  taking the solution to be determined by
C
C                        2     2
C              {\phi} = x  -  y        (in the x-y plane).
C
C The form of the test problems are illustrated in the following
C  diagram.
C
C         ^
C       y |
C         |                    ^n
C                              |
C        0.1   --------------------------------------
C              |B                                  C|
C              |                                    |
C              |       *                     *      |
C              |               2     2              |
C          <-  |     {\phi} = x  -  y               | ->
C           n  |                                    |  n
C              |                  *                 |
C              |                                    |
C              |                                    |         
C              |                                    |
C              |        *                   *       |
C              |                                    |
C              |A                                  D|
C        0.0   --------------------------------------
C                              | n
C              0.0            \ /                0.1  ---> x
C
C Note that on side CD  v = d {\phi}/dx, on AB  v = - d {\phi}/dx,
C  on BC  v = d {\phi}/dy and on DA  v = - d{\phi}/dy: the sign
C  is negative if the normal is in the negative x or y direction.
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 interior points at which the solution is sought;
C  the points (0.025,0.025), (0.075,0.025), (0.025,0.075),
C  (0.075,0.075), and (0.05,0.05).

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 MAXNPI  : The limit on the number of interior points.


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


C The program 

      PROGRAM LIBEM2T
      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 interior to the boundary, where 
C    properties are sought
      INTEGER    MAXNPI
      PARAMETER (MAXNPI=6)

C  Constants
C   Real scalars: 0, 1, 2, pi
      REAL*8 ZERO,ONE,TWO,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 interior to the boundary(ies) where the 
C    properties are sought and the directional vectors at those points.
C    [Only necessary if an interior solution is sought.]
C    Number of interior points and counter
      INTEGER    NPI,IPI
C    Coordinates of the interior points
      REAL*8     PINT(MAXNPI,2)


C   Data structures that are used to define each test problem in turn
C    and are input parameters to LIBEM2.
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    SIPHI(j) is the value of the free-field potential at the centre 
C     of the j-th panel.
      REAL*8 SIPHI(MAXNS)
C    PDIPHI(j) is the value of the free-field potential at the interior
C     points.
      REAL*8 PDIPHI(MAXNPI)


C  Validation and control parameters for LIBEM2
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 LIBEM2
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 interior points
      REAL*8 PDPHI(MAXNPI)

C Workspace for LIBEM2
      REAL*8  L_SS(MAXNS,MAXNS)
      REAL*8  M_SS(MAXNS,MAXNS)
      REAL*8  L_DS(MAXNPI,MAXNS)
      REAL*8  M_DS(MAXNPI,MAXNS)
      REAL*8  C(MAXNS)
      REAL*8  WKSPC1(MAXNS)
      REAL*8  WKSPC2(MAXNS)
      REAL*8  WKSPC3(MAXNS)
      INTEGER*4 PERM(MAXNS)
      LOGICAL XORY(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  The normals to the panels  
      REAL*8     NORMAL(MAXNS,2)
C  Functions defining the exact {\phi} and d{\phi}/dn
      REAL*8     PHI,DPHIDN
C  Temporary variables
      REAL*8     Q(2),QNORM(2),QA(2),QB(2)

C  Other variables
      REAL*8     EPS


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

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      PI=4.0D0*ATAN(ONE)
      EPS=1.0E-10



C Describe the nodes and panels that make up the boundary
C  :The square with vertices (0,0), (0.1,0), (0.1,0.1), (0,0.1) is divided 
C  : into NS=32 uniform panels. VERTEX and SELV are defined 
C  : anti-clockwise around the boundary so that the normal to the 
C  : boundary is assumed to be outward
C  :Set up nodes
C  : Set NS, the number of panels
      NS=32
C  : Set NV, the number of vertices (equal to the number of panels)
      NV=NS
C  : Set coordinates of the nodes
      DATA ((VERTEX(IV,ICOORD),ICOORD=1,2),IV=1,32)
     * / 0.0000000000D0, 0.00000000000D0,
     *   0.0000000000D0, 0.1250000000D0,
     *   0.0000000000D0, 0.2500000000D0,
     *   0.0000000000D0, 0.3750000000D0,
     *   0.0000000000D0, 0.5000000000D0,
     *   0.0000000000D0, 0.6250000000D0,
     *   0.0000000000D0, 0.7500000000D0,
     *   0.0000000000D0, 0.8750000000D0,

     *   0.0000000000D0, 1.0000000000D0,
     *   0.1250000000D0, 1.0000000000D0,
     *   0.2500000000D0, 1.0000000000D0,
     *   0.3750000000D0, 1.0000000000D0,
     *   0.5000000000D0, 1.0000000000D0,
     *   0.6250000000D0, 1.0000000000D0,
     *   0.7500000000D0, 1.0000000000D0,
     *   0.8750000000D0, 1.0000000000D0,

     *   1.0000000000D0, 1.0000000000D0,
     *   1.0000000000D0, 0.8750000000D0,
     *   1.0000000000D0, 0.7500000000D0,
     *   1.0000000000D0, 0.6250000000D0,
     *   1.0000000000D0, 0.5000000000D0,
     *   1.0000000000D0, 0.3750000000D0,
     *   1.0000000000D0, 0.2500000000D0,
     *   1.0000000000D0, 0.1250000000D0,

     *   1.00000000000D0, 0.0000000000D0,
     *   0.8750000000D0, 0.0000000000D0,
     *   0.7500000000D0, 0.0000000000D0,
     *   0.6250000000D0, 0.0000000000D0,
     *   0.5000000000D0, 0.0000000000D0,
     *   0.3750000000D0, 0.0000000000D0,
     *   0.2500000000D0, 0.0000000000D0,
     *   0.1250000000D0, 0.0000000000D0 /

C  :Describe the panels that make up the two boundarys
C  : Set NS, the number of panels
      NS=32
C  : Set nodal indices that describe the panels of the boundary(s).
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,32)
     * /  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,
     *   19,  20,
     *   20,  21,
     *   21,  22,
     *   22,  23,
     *   23,  24,
     *   24,  25,

     *   25,  26,
     *   26,  27,
     *   27,  28,
     *   28,  29,
     *   29,  30,
     *   30,  31,
     *   31,  32,
     *   32,  1 /
       
C Set the centres of the panels, the central points
      DO 100 IS=1,NS
        QA(1)=VERTEX(SELV(IS,1),1)
        QA(2)=VERTEX(SELV(IS,1),2)
        QB(1)=VERTEX(SELV(IS,2),1)
        QB(2)=VERTEX(SELV(IS,2),2)
        SELCNT(IS,1)=(QA(1)+QB(1))/TWO
        SELCNT(IS,2)=(QA(2)+QB(2))/TWO
        CALL NORM2(QA,QB,QNORM)
        NORMAL(IS,1)=QNORM(1)
        NORMAL(IS,2)=QNORM(2)
100   CONTINUE


C Set the points in the domain where the  properties
C  are sought, PINT. 
C : Let them be the points (0.025,0.025), (0.075,0.025), (0.025,0.075),
C :  (0.075,0.075) and (0.5,0.5).
      NPI=5
      DATA ((PINT(IPI,ICOORD),ICOORD=1,2),IPI=1,5)
     *  /  0.250D0,     0.250D0,
     *     0.750D0,     0.250D0,
     *     0.250D0,     0.750D0,
     *     0.750D0,     0.750D0,
     *     0.500D0,     0.500D0 /


C The number of points on the boundary is equal to the number of 
C  panels
      NSP=NS
        
C Set up test problems

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

C   Set up particular alpha and beta functions for the boundary 
C    condition
C  :Set nature of the boundary condition by prescribing the values of
C   the boundary functions SALVAL and SBEVAL at the central points
      DO 640 ISP=1,NSP/2
C   :In this case a Dirichlet (phi-valued) boundary condition
        SALPHA(ISP)=ONE
        SBETA(ISP)=ZERO
        Q(1)=SELCNT(ISP,1)
        Q(2)=SELCNT(ISP,2)
        SF(ISP)=PHI(Q)
        SIPHI(ISP)=ZERO
640   CONTINUE
      DO 645 ISP=NSP/2+1,NSP
C   :In this case a Neumann boundary condition
        SALPHA(ISP)=ZERO
        SBETA(ISP)=ONE
        Q(1)=SELCNT(ISP,1)
        Q(2)=SELCNT(ISP,2)
        QNORM(1)=NORMAL(ISP,1)
        QNORM(2)=NORMAL(ISP,2)
        SF(ISP)=DPHIDN(Q,QNORM)
        SIPHI(ISP)=ZERO
        C(ISP)=SIPHI(ISP)
645   CONTINUE
      DO 650 IPI=1,NPI
        PDIPHI(IPI)=ZERO
650   CONTINUE

      CALL LIBEM2(MAXNV,NV,VERTEX,MAXNS,NS,SELV,
     *                 MAXNPI,NPI,PINT,
     *                 SALPHA,SBETA,SF,SIPHI,PDIPHI,
     *                 LSOL,LVALID,EGEOM,
     *                 SPHI,SVEL,PDPHI,
     *                 L_SS,M_SS,L_DS,M_DS,
     *                 C,PERM,XORY,WKSPC1,WKSPC2,WKSPC3)
      

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

C  Formats for output

      WRITE(20,*) 'LIBEM2: Solution to the test problem LIBEM2_T'
      WRITE(20,*)
C  Loop(ISP) through the points on the boundary
      WRITE(20,*) '     x         y          phi_app   phi_ex'
      DO 210 IPI=1,NPI
        Q(1)=PINT(IPI,1)
        Q(2)=PINT(IPI,2)
        WRITE(20,999)  Q(1),Q(2),PDPHI(IPI),PHI(Q)
210   CONTINUE  

C   Close file
      CLOSE(20)

999   FORMAT(2F10.3,'    ',2F10.5)
     

      END


C  PHI = x**2 - y**2
      REAL*8 FUNCTION PHI(Q)
      REAL*8 Q(2)
      PHI=2*(Q(1)**2-Q(2)**2)
      END

      REAL*8 FUNCTION DPHIDN(Q,QNORM)
      REAL*8 Q(2),QNORM(2)
      DPHIDN=2*(2*Q(1)*QNORM(1)-2*Q(2)*QNORM(2))
      END
                                                                          