C*************************************************************     
C Test program for subroutine LIBEM3 by Stephen Kirkup       *    
C************************************************************* 
C
C  Version 2 (Version 1 2001). Revised test problems. Changes
C   in LIBEM3.
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/LIBEM3_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 LIBEM3. The program computes
C  the solution to a Laplace problem interior to a sphere by the 
C  boundary panel method.
C
C Background
C ----------
C
C We wish to solve the Laplace equation
C
C                  __ 2              
C                  \/    {\phi}   =  0  
C
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 real-valued functions defined
C   on S. 
C
C Subroutine LIBEM3 accepts a description of the boundary of the domain 
C  and the position of the interior 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 interior points.
C

C The test problems
C -----------------
C
C In this test the domain is a sphere of radius 1 (metre). The solution
C  to the problem with a Dirichlet boundary condition ({\alpha}=1, 
C  {\beta}=0) on the upper surface (z>=0) and with a Neumann boundary 
C  condition  ({\alpha}=0, beta=1) on the lower surface (z<0) are 
C  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} = x   +   y   +   z  ,
C
C  which is clearly a solution of the Laplace equation.
C
C The boundary is described by a set of NS=36 planar triangular panels
C  of approximately equal size. The boundary solution points are the 
C  centres of the panels. 
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 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 LIBEM3: Subroutine for solving the interior Laplace
C  equation (file LIBEM3.FOR contains the subroutine LIBEM3)
C subroutine H3LC: Returns the individual discrete Laplace integral
C  operators. (file H3LC.FOR contains H3LC and subordinate routines)
C subroutine GLS: Solves a general linear system of equations.
C  (file GLS.FOR contains CGSL and subordinate routines)
C file GEOM3D.FOR contains the set of relevant geometric subroutines


C The program 

      PROGRAM LIBEM3T
      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=40)
C   Limit on the number of vertices (equal to the number of panels)
      INTEGER    MAXNV
      PARAMETER (MAXNV=30)
C   Limit on the number of points interior to the boundary, where 
C    solutions are sought
      INTEGER    MAXNPI
      PARAMETER (MAXNPI=6)

C  Constants
C   Real scalars: 0, 1, 2, pi
      REAL*8 ZERO,ONE,THREE

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,3)
C   The three nodes that define each panel on the boundaries
      INTEGER    SELV(MAXNS,3)
C   The points interior to the boundary(ies) where the solution
C    is 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   Data structures that are used to define each test problem in turn
C    and are input parameters to LIBEM3.
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   Incident field
C    real SIPHI(MAXNS): The incident velocity potential at the
C     centres of the panels
      REAL*8 SIPHI(MAXNS)

C    real PDIPHI(MAXNPI): The incident velocity potential at the chosen
C    interior points
      REAL*8 PDIPHI(MAXNPI)

C  Validation and control parameters for LIBEM3
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 LIBEM3
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  Discrete Operators (L_SS and M_SSPLUS changed in GLS)
      REAL*8     L_SS(MAXNS,MAXNS)
      REAL*8     M_SSPLUS(MAXNS,MAXNS)
      REAL*8     L_DS(MAXNPI,MAXNS)
      REAL*8     M_DS(MAXNPI,MAXNS)

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

C Work Space
      REAL*8     WKSPC1(MAXNS)
      REAL*8     WKSPC2(MAXNS)
      REAL*8     WKSPC3(MAXNS)
      
      REAL*8 PHI,DPHI

      REAL*8 SIZE3
      REAL*8 P(3),NP(3),SIZEP

C  Counter through the x,y coordinates
      INTEGER    ICOORD

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



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

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      THREE=3.0D0




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


C   Set up VERTEX, the coordinates of the vertices of the panels
      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 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,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 panels, the central 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 are sought, PINT. 
C : Set four points
      NPI=4
      DATA ((PINT(IPI,ICOORD),ICOORD=1,3),IPI=1,4)
     *  /  0.500D0,     0.000D0,  0.0000D0,
     *     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  panels
      NSP=NS
        

C  :Switch for particular solution
      LSOL=.TRUE.
C  :Switch on the validation of LIBEM3
      LVALID=.TRUE.
C  :Set EGEOM
      EGEOM=1.0D-6

C  Open file for the output data
      OPEN(UNIT=20,FILE='LIBEM3.OUT')
      WRITE(20,*) 'Tests on LIBEM3'
      WRITE(20,*) 'A subroutine for solving the 3D interior 
     * Laplace Equation'
      WRITE(20,*) 'The surface is the unit sphere that is divided
     * into 36 panels'

C TEST 1

      WRITE(20,*)
      WRITE(20,*) 'Test 1 - Dirichlet - {\phi}=1' 
      WRITE(20,*)
C   Set up particular alpha and beta functions for Dirichlet test 
C    and type of boundary condition
      DO 140 ISP=1,NSP
        SALPHA(ISP)=1.0D0
        SBETA(ISP)=0.0D0
        SF(ISP)=1.0D0
140   CONTINUE
      DO 160 ISP=1,NSP
        SIPHI(ISP)=ZERO
160   CONTINUE
      DO 170 IPI=1,NPI
        PDIPHI(IPI)=ZERO
170   CONTINUE

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


      WRITE(20,*) 'Solutions at the interior points'
      WRITE(20,*) '   P(1)      P(2)      P(3)     Computed    Exact'
      DO 180 IPI=1,NPI
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)
        P(3)=PINT(IPI,3)
        WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),1.0
180   CONTINUE

C TEST 2

      WRITE(20,*)
      WRITE(20,*) 'Test 2 - Half Dirichlet, half Neumann - 
     * {\phi}=1, v=0' 
      WRITE(20,*)
C   Set up particular alpha and beta functions for Dirichlet/Neumann test 
C    and type of boundary condition
      DO 240 ISP=1,NSP/2
        SALPHA(ISP)=1.0D0
        SBETA(ISP)=0.0D0
        SF(ISP)=1.0D0
240   CONTINUE
      DO 250 ISP=NSP/2+1,NSP
        SALPHA(ISP)=0.0D0
        SBETA(ISP)=1.0D0
        SF(ISP)=0.0D0
250   CONTINUE
      DO 260 ISP=1,NSP
        SIPHI(ISP)=ZERO
260   CONTINUE
      DO 270 IPI=1,NPI
        PDIPHI(IPI)=ZERO
270   CONTINUE

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


      WRITE(20,*) 'Solutions at the interior points'
      WRITE(20,*) '   P(1)      P(2)      P(3)     Computed    Exact'
      DO 280 IPI=1,NPI
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)
        P(3)=PINT(IPI,3)
        WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),1.0
280   CONTINUE

C TEST 3

      WRITE(20,*)
      WRITE(20,*) 'Test 3 - Dirichlet - {\phi}=X^2-Y^2+Z' 
      WRITE(20,*)
C   Set up particular alpha and beta functions for Dirichlet/Neumann test 
C    and type of boundary condition
      DO 340 ISP=1,NSP
        P(1)=SELCNT(ISP,1)
        P(2)=SELCNT(ISP,2)
        P(3)=SELCNT(ISP,3)
        SALPHA(ISP)=1.0D0
        SBETA(ISP)=0.0D0
        SF(ISP)=P(1)**2-P(2)**2+P(3)
340   CONTINUE
      DO 360 ISP=1,NSP
        SIPHI(ISP)=ZERO
360   CONTINUE
      DO 370 IPI=1,NPI
        PDIPHI(IPI)=ZERO
370   CONTINUE

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


      WRITE(20,*) 'Solutions at the interior points'
      WRITE(20,*) '   P(1)      P(2)      P(3)     Computed    Exact'
      DO 380 IPI=1,NPI
        P(1)=PINT(IPI,1)
        P(2)=PINT(IPI,2)
        P(3)=PINT(IPI,3)
        WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),P(1)**2-P(2)**2+P(3)
380   CONTINUE



999   FORMAT(5F10.4)
      CLOSE(20)
      END

      REAL*8 FUNCTION PHI(P)
      REAL*8 P(3)
      PHI=P(1)*P(1)-2*P(2)*P(2)+P(3)*P(3)
      END

      REAL*8 FUNCTION DPHI(P,NP)
      REAL*8 P(3),NP(3)
      DPHI=NP(1)+NP(2)+NP(3)
      END



                                                                            