C***************************************************************
C     Test program for subroutine HIBEM3 by Stephen Kirkup     
C***************************************************************
C
C  Copyright 2001- Stephen Kirkup
C  School of Computing Engineering and Physical Sciences
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/HIBEM3_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 HIBEM3. The program computes
C  the solution to Helmholtz problem interior to a sphere by the 
C  boundary element method.
C
C Background
C ----------
C
C The Helmholtz equation is as follows
C                  __ 2                2
C                  \/    {\phi}   +   k  {\phi}   =  0  
C
C  where k is the wavenumber and {\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-valued functions defined
C   on S. 
C
C Subroutine HIBEM3 accepts the wavenumber, 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 sphere of radius 1 (metre). The 
C  solution to the problem with a Dirichlet boundary condition 
C  ({\alpha}=1, {\beta}=0) and with a Neumann boundary condition 
C  ({\alpha}=0, beta=1) are sought.
C Assuming the sphere is centred at the origin the boundary conditions
C  are specified through taking the solution to be determined by
C
C                      {\phi} = sin(k z),
C
C  which is clearly a solution of the Helmholtz equation.
C
C The boundary is described by a set of NS=36 planar triangular elements
C  of approximately equal size. The boundary solution points are the 
C  centres of the elements. 
C The solution is sought at the interior points (0,0), (0,0.5), 
C  (0.25,0.25).

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


C External modules related to the package
C ---------------------------------------
C subroutine HIBEM3: Subroutine for solving the interior Helmholtz
C  equation (file HIBEM3.FOR contains the subroutine HIBEM3)
C subroutine H3LC: Returns the individual discrete Helmholtz integral
C  operators. (file H3LC.FOR contains H3LC and subordinate routines)
C subroutine CGLS: Solves a general linear system of equations.
C  (file CGLS.FOR contains CGSL and subordinate routines)
C file GEOM3D.FOR contains the set of relevant geometric subroutines


C The program 

      PROGRAM IBEM3T
      IMPLICIT NONE

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=40)
C   Limit on the number of vertices (equal to the number of elements)
      INTEGER    MAXNV
      PARAMETER (MAXNV=30)
C   Limit on the number of test problems
      INTEGER    MAXTEST
      PARAMETER (MAXTEST=5)
C   Limit on the number of points interior to the boundary, where 
C    solution values are sought
      INTEGER    MAXNPI
      PARAMETER (MAXNPI=6)

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

C   Wavenumber parameter for HIBEM3
      COMPLEX*16     K

C  Geometrical description of the boundary(ies)
C   Number of elements and counter
      INTEGER    NS,IS
C   Number of collocation 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,3)
C   The three nodes that define each element on the boundaries
      INTEGER    SELV(MAXNS,3)
C   The points interior to the boundary(ies) where the solution 
C    values 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,3)

C  Number of test problems and counter
      INTEGER    NTEST,ITEST

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   Incident field 
C    The  potential at the centres of the elements in each test
C     case
      COMPLEX*16 SPHIIN(MAXTEST,MAXNS)
C    The derivative of the potential at the centres of the 
C     elements in each test
      COMPLEX*16 SVELIN(MAXTEST,MAXNS)
C    The potential at the interior points 
      COMPLEX*16 PPHIIN(MAXTEST,MAXNPI)

C   Storage for the returned results 
C    The  potential at the centres of the elements in each test
C     case
      COMPLEX*16 SPHIS(MAXTEST,MAXNS)
C    The derivative of the potential at the centres of the 
C     elements in each test
      COMPLEX*16 SVELS(MAXTEST,MAXNS)
C    The potential at the interior points 
      COMPLEX*16 PPHIS(MAXTEST,MAXNPI)

C    Storage of exact results
      COMPLEX*16 EXPHI(MAXTEST,MAXNPI)

C   Data structures that are used to define each test problem in turn
C    and are input parameters to HIBEM3.
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(MAXNPI)

C  Validation and control parameters for HIBEM3
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 parameter
      COMPLEX*16 MU

C Output from subroutine HIBEM3
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 interior points
      COMPLEX*16 PIPHI(MAXNPI)

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


C  Counter through the x,y coordinates
      INTEGER    ICOORD

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

C  Other variables used in specifying the boundary condition
      REAL*8     Z,QA(3),QB(3),QC(3),QNORM(3),DZBDN

      REAL*8     EPS


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

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      THREE=3.0D0
      FOUR=4.0D0
      PI=4.0D0*ATAN(ONE)
      CZERO=CMPLX(ZERO,ZERO)
      CONE=CMPLX(ONE,ZERO)
      CIMAG=CMPLX(ZERO,ONE)

      EPS=1.0E-10


C Describe the nodes and elements that make up the boundary
C  :The unit sphere, centred at the origin is divided 
C  : into NS=36 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=36
C  : Set NV, the number of vertices (equal to the number of elements)
      NV=NS
C  : Set coordinates of the nodes


C   Set up VERTEX, the coordinates of the vertices of the elements
      NV=20
      DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,20)
     * / 0.000D0, 0.000D0, 1.000D0,
     *   0.000D0, 0.745D0, 0.667D0,
     *   0.645D0, 0.372D0, 0.667D0,
     *   0.645D0,-0.372D0, 0.667D0,
     *   0.000D0,-0.745D0, 0.667D0,
     *  -0.645D0,-0.372D0, 0.667D0,
     *  -0.645D0, 0.372D0, 0.667D0,
     *   0.500D0, 0.866D0, 0.000D0,
     *   1.000D0, 0.000D0, 0.000D0,
     *   0.500D0,-0.866D0, 0.000D0,
     *  -0.500D0,-0.866D0, 0.000D0,
     *  -1.000D0, 0.000D0, 0.000D0,
     *  -0.500D0, 0.866D0, 0.000D0,
     *   0.000D0, 0.745D0,-0.667D0,
     *   0.645D0, 0.372D0,-0.667D0,
     *   0.645D0,-0.372D0,-0.667D0,
     *   0.000D0,-0.745D0,-0.667D0,
     *  -0.645D0,-0.372D0,-0.667D0,
     *  -0.645D0, 0.372D0,-0.667D0,
     *   0.000D0, 0.000D0,-1.000D0/


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 acoustic domain.
      DATA ((SELV(IS,ICOORD),ICOORD=1,3),IS=1,36)
     * / 1, 3, 2,    1, 4, 3,    1, 5, 4,    1, 6, 5,
     *   1, 7, 6,    1, 2, 7,    2, 3, 8,    3, 9, 8,
     *   3, 4, 9,    4,10, 9,    4, 5,10,    5,11,10,
     *   5, 6,11,    6,12,11,    6, 7,12,    7,13,12,
     *   7, 2,13,    2, 8,13,    8,15,14,    8, 9,15,
     *   9,16,15,    9,10,16,   10,17,16,   10,11,17,
     *  11,18,17,   11,12,18,   12,19,18,   12,13,19,
     *  13,14,19,   13, 8,14,   14,15,20,   15,16,20,
     *  16,17,20,   17,18,20,   18,19,20,   19,14,20/


C Set the centres of the elements, the collocation points
      DO IS=1,NS
        SELCNT(IS,1)=(VERTEX(SELV(IS,1),1)
     *   +VERTEX(SELV(IS,2),1)+VERTEX(SELV(IS,3),1))/THREE
        SELCNT(IS,2)=(VERTEX(SELV(IS,1),2)
     *   +VERTEX(SELV(IS,2),2)+VERTEX(SELV(IS,3),2))/THREE
        SELCNT(IS,3)=(VERTEX(SELV(IS,1),3)
     *   +VERTEX(SELV(IS,2),3)+VERTEX(SELV(IS,3),3))/THREE
      END DO



C Set the points in the domain where the solutions
C  are sought, PINT. 
C : Let them be the points (0.000,0.000), (0.000,0.500), (0.250,0.250).
      NPI=5
      DATA ((PINT(IPI,ICOORD),ICOORD=1,3),IPI=1,5)
     *  /  0.500D0,     0.000D0,  0.0000D0,
     *     0.000D0,     0.000D0,  0.010D0,
     *     0.000D0,     0.000D0,  0.250D0,
     *     0.000D0,     0.000D0,  0.500D0,
     *     0.000D0,     0.000D0,  0.750D0/


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=2


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

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

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 
C    {\phi}=sin((k z)
C   :Set K, the wavenumber
      K=KVAL(1)
      DO 170 ISP=1,NSP
        Z=SELCNT(ISP,3)
        SFVAL(1,ISP)=SIN(K*Z)
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,NPI
        PPHIIN(1,IPI)=0.0D0
190   CONTINUE

C The exact solution at the interior points
      DO 195 IPI=1,NPI
        Z=PINT(IPI,3)
        EXPHI(1,IPI)=SIN(K*Z)
195   CONTINUE



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

C  : Set the wavenumber in KVAL
      KVAL(2)=1.5D0

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 200 ISP=1,NSP
        SALVAL(2,ISP)=CZERO
        SBEVAL(2,ISP)=CONE
200   CONTINUE

C  :The test problem is devised so that 
C    {\phi}=sin((k z)
C   :Set K, the wavenumber
      K=KVAL(2)
      DO 210 ISP=1,NSP
        Z=SELCNT(ISP,3)
        QA(1)=VERTEX(SELV(ISP,1),1)
        QA(2)=VERTEX(SELV(ISP,1),2)
        QA(3)=VERTEX(SELV(ISP,1),3)
        QB(1)=VERTEX(SELV(ISP,2),1)
        QB(2)=VERTEX(SELV(ISP,2),2)
        QB(3)=VERTEX(SELV(ISP,2),3)
        QC(1)=VERTEX(SELV(ISP,3),1)
        QC(2)=VERTEX(SELV(ISP,3),2)
        QC(3)=VERTEX(SELV(ISP,3),3)
        CALL NORM3(QA,QB,QC,QNORM)
        DZBDN=QNORM(3)
        SFVAL(2,ISP)=K*COS(K*Z)*DZBDN
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,NPI
        EXPHI(2,IPI)=0.0D0
230   CONTINUE

C The exact solution at the interior points
      DO 295 IPI=1,NPI
        Z=PINT(IPI,3)
        EXPHI(2,IPI)=SIN(K*Z)
295   CONTINUE




C  :Switch for particular solution
      LSOL=.TRUE.
C  :Switch on the validation of HIBEM3
      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 IPI=1,NPI
            PFFPHI(IPI)=PPHIIN(ITEST,IPI)
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 HIBEM3(K,
     *                 MAXNV,NV,VERTEX,MAXNS,NS,SELV,
     *                 MAXNPI,NPI,PINT,
     *                 SALPHA,SBETA,SF,SFFPHI,SFFVEL,PFFPHI,
     *                 LSOL,LVALID,EGEOM,MU,
     *                 SPHI,SVEL,PIPHI,
     *                 WKSPC1,WKSPC2,WKSPC3,WKSPC4,
     *                 WKSPC5,WKSPC6,WKSPC7)

C   Set up particular alpha and beta functions for this wavenumber
C    and type of boundary condition
          DO 660 ISP=1,NSP
            SPHIS(ITEST,ISP)=SPHI(ISP)
            SVELS(ITEST,ISP)=SVEL(ISP)
660       CONTINUE
          DO 670 IPI=1,NPI
            PPHIS(ITEST,IPI)=PiPHI(IPI)
670       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='HIBEM3.OUT')

C  Formats for output
2810  FORMAT(1X,'Wavenumber = ',F8.2/)
2845  FORMAT(4X,'Potential at the interior points',/)
2850  FORMAT(5X,'index',7X,'Potential',19X,'Exact Potential'/)
2860  FORMAT(4X,I4,2X,E10.4,'+ ',E10.4,' i    ',
     * E10.4, '+ ',E10.4,' i    '/)

      WRITE(20,*) 'HIBEM3: Computed solution to the Helmholtz problems'
      WRITE(20,*)
C  Loop(ITEST) through the test problems.
      DO 2000 ITEST=1,NTEST
C   Output the wavenmber
        WRITE(20,*)
        WRITE(20,*)
        WRITE(20,*) 'Test problem ',ITEST
        WRITE(20,*)
        WRITE(20,2810) KVAL(ITEST)
        WRITE(20,*)
        WRITE(20,2845)
        WRITE(20,*)
        WRITE(20,2850)
        DO 2020 IPI=1,NPI
        WRITE(20,2860) IPI,DREAL(PPHIS(ITEST,IPI)), 
     *  DIMAG(PPHIS(ITEST,IPI)), 
     *  DREAL(EXPHI(ITEST,IPI)), 
     *  DIMAG(EXPHI(ITEST,IPI))
2020    CONTINUE
        WRITE(20,*)
        WRITE(20,*)
        WRITE(20,*)
2000  CONTINUE  

      CLOSE(20)

      END



                                                                                       