C***************************************************************
C     Test program for subroutine HEBEM2 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/HEBEM2_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 HEBEM2. The program computes
C  the solution to an acoustic/Helmholtz problem exterior to a square 
C  by the boundary element method.
C

C In this test the domain is a exterior to a square of side 0.1 (metre).
C Assuming the boundary condition is derived through assuming it is 
C  equivalent to a point source at the centre of the square is clearly 
C  a solution of the Helmholtz equation.
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              |                                    |
C          <-  |                                    | ->
C           n  |                                    |  n
C              |                  *                 |
C              |           equivalent source        |
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 elements of equal size,
C  so that each side comprises eight elements. The boundary solution
C  points are the centres of the elements. 
C The * shows the position of the equivalent source.

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 elements.
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 in the package
C ---------------------------------------
C subroutine HEBEM2: Subroutine for solving the exterior Helmholtz
C  equation (file HEBEM2.FOR contains IBEM2 and subordinate routines)
C subroutine H2LC: Returns the individual discrete Helmholtz integral
C  operators. (file H2LC.FOR contains H2LC and subordinate routines)
C subroutine CGLS: Solves a general linear system of equations.
C  (file CGLS.FOR contains CGSL and subordinate routines)
C subroutine FNHANK: This computes Hankel functions of the first kind
C  and of order zero and one. (e.g. file FNHANK.FOR)
C file GEOM2D contains the set of relevant geometric subroutines

C The program 

      PROGRAM  HEBEM2T
      IMPLICIT REAL*8 (A-H,O-Z)

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

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

C  Constants
C   Real scalars: 0, 1, 2, pi
      REAL*8 ZERO,ONE,TWO,PI
C   Complex scalars: (0,0), (1,0), (0,1)
      COMPLEX*16 CZERO,CONE,CIMAG

C   Wavenumber parameter for HEBEM2
      COMPLEX*16     K

C  Geometrical description of the boundary(ies)
C   Number of elements and counter
      INTEGER    NS,IS
C   Number of collocation points 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 element on the boundaries
      INTEGER    SELV(MAXNS,2)
C   The points exterior to the boundary(ies) where the 
C    properties are sought and the directional vectors at those points.
C    [Only necessary if an exteriorsolution is sought.]
C    Number of exterior points
      INTEGER    NPE
C    Coordinates of the exterior points
      REAL*8     PEXT(MAXNPE,2)

C  Data structures that contain the parameters that define the test
C   problems
C   The wavenumber for each test. KVAL(i) is assigned the wavenumber
C    of the i-th test problem.
      REAL*8     KVAL(MAXTEST)
C   The nature of the boundary condition is specified by assigning 
C    values to the data structures SALVAL and SBEVAL. 
C    SALVAL(i,j) is assigned the value of {\alpha} at the center of the
C     j-th element for the i-th test problem.
      COMPLEX*16 SALVAL(MAXTEST,MAXNS)
C    SBEVAL(i,j) is assigned the value of {\beta} at the center of the
C     j-th element for the i-th test problem.
      COMPLEX*16 SBEVAL(MAXTEST,MAXNS)      
C   The actual boundary condition is specified by assigning values to 
C    the data structure SFVAL. 
C    SFVAL(i,j) is assigned the value of f at the center of the j-th 
C    element for the i-th test problem.
      COMPLEX*16 SFVAL(MAXTEST,MAXNS)
C    The incident potential at the centres of the elements in each test
C     case
      COMPLEX*16 SPHIIN(MAXTEST,MAXNS)
C    The derivative of the incident potential at the centres of the 
C     elements in each test
      COMPLEX*16 SVELIN(MAXTEST,MAXNS)
C    The incident potential at the exterior points 
      COMPLEX*16 PPHIIN(MAXTEST,MAXNPE)

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

C    The incident potential at the centres of the elements in each test
C     case
      COMPLEX*16 SFFPHI(MAXNS)
C    The derivative of the incident potential at the centres of the 
C     elements in each test
      COMPLEX*16 SFFVEL(MAXNS)
C    The incident potential at the interior points 
      COMPLEX*16 PFFPHI(MAXNPE)

C  Validation and control parameters for HEBEM2
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   The weighting parameter.
      COMPLEX*16 MU

C Output from subroutine HEBEM2
C  The velocity potential (phi - the solution) at the centres of the 
C   elements
      COMPLEX*16 SPHI(MAXNS)
C  The normal derivative of the velocity potential at the centres of the
C    elements
      COMPLEX*16 SVEL(MAXNS)
C  The velocity potential (phi - the solution) at exterior points
      COMPLEX*16 PEPHI(MAXNPE)

C Workspace for HEBEM2
      COMPLEX*16 WKSPC1(MAXNS,MAXNS)
      COMPLEX*16 WKSPC2(MAXNS,MAXNS)
      COMPLEX*16 WKSPC3(MAXNPE,MAXNS)
      COMPLEX*16 WKSPC4(MAXNPE,MAXNS)
      COMPLEX*16 WKSPC5(MAXNS)
      COMPLEX*16 WKSPC6(MAXNS)
      LOGICAL    WKSPC7(MAXNS)

C Number of test problems and counter
      INTEGER    NTEST,ITEST


C   Helmholtz properties. These data structures are appended after each
C    execution of HEBEM2 and contain the numerical solution to the test
C    problems. 
C    At the centres of the elements
C     Velocity potential phi on the boundary
      COMPLEX*16 SPHIVAL(MAXTEST,MAXNS)
C     Velocity potential phi at the exterior points
      COMPLEX*16 EPHIVAL(MAXTEST,MAXNPE)

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

C  Variables used in setting up test boundary conditions
      REAL*8     CENTX,CENTY,DRBDN,RR(2)
      COMPLEX*16 H(0:1),KR

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)
      CZERO=CMPLX(ZERO,ZERO)
      CONE=CMPLX(ONE,ZERO)
      CIMAG=CMPLX(ZERO,ONE)

      EPS=1.0E-10

C  Reference for decibel scales
      PREREF=2.0D-05


C Describe the nodes and elements 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 elements. 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 elements
      NS=32
C  : Set NV, the number of vertices (equal to the number of elements)
      NV=NS
C  : Set coordinates of the nodes
      DATA ((VERTEX(IV,I),I=1,2),IV=1,32)
     * / 0.00000000000D0, 0.00000000000D0,
     *   0.00000000000D0, 0.01250000000D0,
     *   0.00000000000D0, 0.02500000000D0,
     *   0.00000000000D0, 0.03750000000D0,
     *   0.00000000000D0, 0.05000000000D0,
     *   0.00000000000D0, 0.06250000000D0,
     *   0.00000000000D0, 0.07500000000D0,
     *   0.00000000000D0, 0.08750000000D0,

     *   0.00000000000D0, 0.10000000000D0,
     *   0.01250000000D0, 0.10000000000D0,
     *   0.02500000000D0, 0.10000000000D0,
     *   0.03750000000D0, 0.10000000000D0,
     *   0.05000000000D0, 0.10000000000D0,
     *   0.06250000000D0, 0.10000000000D0,
     *   0.07500000000D0, 0.10000000000D0,
     *   0.08750000000D0, 0.10000000000D0,

     *   0.10000000000D0, 0.10000000000D0,
     *   0.10000000000D0, 0.08750000000D0,
     *   0.10000000000D0, 0.07500000000D0,
     *   0.10000000000D0, 0.06250000000D0,
     *   0.10000000000D0, 0.05000000000D0,
     *   0.10000000000D0, 0.03750000000D0,
     *   0.10000000000D0, 0.02500000000D0,
     *   0.10000000000D0, 0.01250000000D0,

     *   0.10000000000D0, 0.00000000000D0,
     *   0.08750000000D0, 0.00000000000D0,
     *   0.07500000000D0, 0.00000000000D0,
     *   0.06250000000D0, 0.00000000000D0,
     *   0.05000000000D0, 0.00000000000D0,
     *   0.03750000000D0, 0.00000000000D0,
     *   0.02500000000D0, 0.00000000000D0,
     *   0.01250000000D0, 0.00000000000D0 /

C  :Describe the elements that make up the two boundaries
C  : Set NS, the number of elements
      NS=32
C  : Set nodal indices that describe the elements of the boundaries.
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(IE,J),J=1,2),IE=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 elements, the collocation points
      DO 100 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
100   CONTINUE


C Set the points in the acoustic domain where the properties
C  are sought, PEXT. 
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).
      NPE=4
      DATA ((PEXT(IPE,ICOORD),ICOORD=1,2),IPE=1,4)
     *  /  0.0000D0,     0.1500D0,
     *     0.0500D0,     0.1500D0,
     *     0.1000D0,     0.1500D0,
     *     0.0500D0,    -0.1000D0 /


C The number of points on the boundary is equal to the number of 
C  elements
      NSP=NS
        
C Set up test problems
C  :Set the number of test problems
      NTEST=3


C  TEST PROBLEM 1
C  ==============

C  : Set the wavenumber in KVAL
      KVAL(1)=CONE

C  :Set nature of the boundary condition by prescribing the values of
C   the boundary functions SALVAL and SBEVAL at the collocation points
C   :In this case a Dirichlet (phi-valued) boundary condition
      DO 160 ISP=1,NSP
        SALVAL(1,ISP)=CONE
        SBEVAL(1,ISP)=CZERO
160   CONTINUE

C  :The test problem is devised so that \phi=i*h0(kr)/4 
C    where r is the distance between the collocation point
C    and the point at the centre of the square
C   :Set K, the wavenumber
      K=KVAL(1)
C   :Set centre of the square
      CENTX=0.05D0
      CENTY=0.05D0
      DO 170 ISP=1,NSP
        X=SELCNT(ISP,1)
        Y=SELCNT(ISP,2)
        R=SQRT((X-CENTX)**2+(Y-CENTY)**2)
        KR=K*R
        CALL FNHANK(KR,H)
        SFVAL(1,ISP)=CIMAG*H(0)/4.0D0
170   CONTINUE     

C  :There is no incident field 
      DO 180 ISP=1,NSP
        SPHIIN(1,ISP)=0.0D0
        SVELIN(1,ISP)=0.0D0
180   CONTINUE
      DO 190 IPI=1,NPE
        PPHIIN(1,IPI)=0.0D0
190   CONTINUE


C  TEST PROBLEM 2
C  ==============

C  : Set the wavenumber in KVAL
      KVAL(2)=CONE

C  :Set nature of the boundary condition by prescribing the values of
C   the boundary functions SALVAL and SBEVAL at the collocation points
C   :In this case a Neumann (v-valued) boundary condition
      DO 200 ISP=1,NSP
        SALVAL(2,ISP)=CZERO
        SBEVAL(2,ISP)=CONE
200   CONTINUE

C   :Set K, the wavenumber
      K=KVAL(2)
C   :Set centre of the square
      CENTX=0.05D0
      CENTY=0.05D0
      DO 210 ISP=1,NSP
        X=SELCNT(ISP,1)
        Y=SELCNT(ISP,2)
        RR(1)=X-CENTX
        RR(2)=Y-CENTY
        R=SQRT(RR(1)**2+RR(2)**2)
        KR=K*R
        CALL FNHANK(KR,H)
C    Face X=0. Note the normal is in the direction of -x
        IF (X.LT.EPS) THEN
         DRBDN=-RR(1)/R
         SFVAL(2,ISP)=-CIMAG*K*DRBDN*H(1)/4.0D0
C    Face X=1. Note that the normal is in the direction of +x
        ELSE IF (X.GT.0.1D0-EPS) THEN
         DRBDN=RR(1)/R
         SFVAL(2,ISP)=-CIMAG*K*DRBDN*H(1)/4.0D0
C    Face Y=0. Note the normal is in the direction of -y
        ELSE IF (Y.LT.EPS) THEN
         DRBDN=-RR(2)/R
         SFVAL(2,ISP)=-CIMAG*K*DRBDN*H(1)/4.0D0
        ELSE 
C    Face Y=1. Note the normal is in the direction of +y
         DRBDN=RR(2)/R
         SFVAL(2,ISP)=-CIMAG*K*DRBDN*H(1)/4.0D0
        END IF
210   CONTINUE     

C   :There is no incident field
      DO 220 ISP=1,NSP
        SPHIIN(2,ISP)=0.0D0
        SVELIN(2,ISP)=0.0D0
220   CONTINUE
      DO 230 IPI=1,NPE
        PPHIIN(2,IPI)=0.0D0
230   CONTINUE


C  TEST PROBLEM 3
C  ==============

C  : Set the wavenumber in KVAL
      KVAL(3)=CONE

C  :The test problem computes the field produced by a unit source at
C    the point (0.5,0.25) within the square with a rigid boundary.
C  :Set nature of the boundary condition by prescribing the values of
C   the boundary functions SALVAL and SBEVAL at the collocation points
C   :The rigid boundary implies the bondary condition v=0.
      DO 240 ISP=1,NSP
        SALVAL(3,ISP)=CZERO
        SBEVAL(3,ISP)=CONE
        SFVAL(3,ISP)=CZERO
240   CONTINUE

C   The incident velocity potential is given by {\phi}_inc=i*h0(kr)/4
C    where r is the distance from the point (0.5,0.25)
      PX=0.05D0
      PY=-0.05D0
C   :Set K, the wavenumber
      K=KVAL(3)
      DO 260 ISP=1,NSP
        X=SELCNT(ISP,1)
        Y=SELCNT(ISP,2)
        RR(1)=X-PX
        RR(2)=Y-PY
        R=SIZE2(RR)
        KR=K*R
        CALL FNHANK(KR,H)
        SPHIIN(3,ISP)=CIMAG*H(0)/4.0D0
C    Face X=0. Note the normal is in the direction of -x
        IF (X.LT.EPS) THEN
          SVELIN(3,ISP)=(-CIMAG*K*H(1)/4.0D0)*(-RR(1)/R)
C    Face X=1. Note that the normal is in the direction of +x
        ELSE IF (X.GT.0.1D0-EPS) THEN
          SVELIN(3,ISP)=(-CIMAG*K*H(1)/4.0D0)*(RR(1)/R)
C    Face Y=0. Note the normal is in the direction of -y
        ELSE IF (Y.LT.EPS) THEN
          SVELIN(3,ISP)=(-CIMAG*K*H(1)/4.0D0)*(-RR(2)/R)
        ELSE 
C    Face Y=1. Note the normal is in the direction of +y
          SVELIN(3,ISP)=(-CIMAG*K*H(1)/4.0D0)*(RR(2)/R)
        END IF
260   CONTINUE     
      DO 270 IPI=1,NPE
        X=PEXT(IPI,1)
        Y=PEXT(IPI,2)
        RR(1)=X-PX
        RR(2)=Y-PY
        R=SIZE2(RR)
        KR=K*R
        CALL FNHANK(KR,H)
        PPHIIN(3,IPI)=CIMAG*H(0)/4.0D0
270   CONTINUE



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

C Loop(ITEST) through the test problems
      DO 500 ITEST=1,NTEST
C  Set OMEGA, the angular frequency omega and K, the wavenumber
        K=KVAL(ITEST)

C   Set up particular alpha and beta functions for this wavenumber
C    and type of boundary condition
          DO 640 ISP=1,NSP
            SALPHA(ISP)=SALVAL(ITEST,ISP)
            SBETA(ISP)=SBEVAL(ITEST,ISP)
            SF(ISP)=SFVAL(ITEST,ISP)
            SFFPHI(ISP)=SPHIIN(ITEST,ISP)
            SFFVEL(ISP)=SVELIN(ITEST,ISP)
640       CONTINUE
          DO 650 IPE=1,NPE
            PEPHI(IPE)=PPHIIN(ITEST,IPE)
650       CONTINUE

C   Set the parameter
          IF(DIMAG(K).GT.DREAL(K)) THEN
            MU=CZERO
          ELSE
           MU=CIMAG/(DREAL(K)-DIMAG(K)+ONE)
          END IF

          CALL HEBEM2(K,
     *                 MAXNV,NV,VERTEX,MAXNS,NS,SELV,
     *                 MAXNPE,NPE,PEXT,
     *                 SALPHA,SBETA,SF,SFFPHI,SFFVEL,PFFPHI,
     *                 LSOL,LVALID,EGEOM,MU,
     *                 SPHI,SVEL,PEPHI,
     *                 WKSPC1,WKSPC2,WKSPC3,WKSPC4,
     *                 WKSPC5,WKSPC6,WKSPC7)


            DO 690 ISP=1,NSP
              SPHIVAL(ITEST,ISP)=SPHI(ISP)
690         CONTINUE

            DO 695 ISP=1,NPE
              EPHIVAL(ITEST,ISP)=PEPHI(ISP)
695         CONTINUE

C  Close loop(ITEST) through the test problems
500     CONTINUE


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

C  Formats for output
2810  FORMAT(1X,'Wavenumber = ',F8.2/)

2860  FORMAT(4X,I4,2X,E10.4,'+ ',E10.4,' i    ')

      WRITE(20,*) 'HEBEM2: Computed solution to the
     * Helmholtz properties'
      WRITE(20,*)
C  Loop(ITEST) through the test problems.
      DO 2000 ITEST=1,NTEST
C   Output the wavenumber
        WRITE(20,*)
        WRITE(20,*)
        WRITE(20,2810) KVAL(ITEST)
        WRITE(20,*)
        WRITE(20,*)
C      Loop(ISP) through the points on the boundary
        DO 2030 ISP=1,NSP
C       Output the velocity potential
          WRITE(20,2860) ISP,DBLE(SPHIVAL(ITEST,ISP)),
     *     AIMAG(SPHIVAL(ITEST,ISP))
2030    CONTINUE
        WRITE(20,*)
C      Loop(IPE) through the points in the exterior
        DO 2040 IPE=1,NPE
C       Output the velocity potential
          WRITE(20,2860) IPE,DBLE(EPHIVAL(ITEST,IPE)),
     *     AIMAG(EPHIVAL(ITEST,IPE))

2040    CONTINUE
        WRITE(20,*)
2000  CONTINUE  

      CLOSE(20)

      END

                                                                                     