C***************************************************************
C     Test program for subroutine LSEM3 by Stephen Kirkup      *        
C***************************************************************
C 
C  Version 2 (July 2015) (Version 1 2001) Changes 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/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 This program is a test for the subroutine LSEM3. The program computes
C  the solution to a Laplace problem exterior to two square plates by 
C  the shell panel method. One plate is assigned the potential one 
C  and the other is at zero potential. The problem may be thought of as 
C  a pair of plates of a capacitor with one plate at zero volts, the 
C  other at 1 volt.
C
C Background
C ----------
C
C We wish to solve the Laplace equation
C
C                  __ 2              
C                  \/    {\phi}   =  0  
C
C
C For the exterior problem, the domain lies exterior to the shells
C  or plates. The shell condition may be Dirichlet, Robin or 
C  Neumann. It is assumed to have the following general form
C
C          a(q) {\delta}(q) + b(q) {\nu}(q) = f(q)
C
C          A(q) {\Phi}(q) + B(q) V(q) = F(q)
C
C  where {\delta}(q)={\phi+}(q) - {\phi-}(q),
C        {\nu}(q) = {\nu+}(q) - {\nu-}(q)
C        {\Phi}(q) = ({\phi+}(q)+{\phi-}(q))/2
C        V(q) = ({v+}(q) + {v-}(q))/2
C  and {\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 LSEM3 accepts a description of the shells and the position
C  of the exterior points where the solution ({\phi}) 
C  is sought, the shell condition and returns the solution ({\delta},
C  {\Phi}, {\nu} and V) on {\Pi} and the value of {\phi} at the exterior 
C  points.
C

C The test problems
C -----------------
C
C In this test the domain is the exterior to two square plates of side
C  0.1m and 0.01m apart. The solution to the problem with a Dirichlet 
C   shell condition (a=1, b=0, f=0) on both plates and (A=1, B=0,
C   F=0) on the lower plate and (A=1, B=0, F=1) on the upper plate.
C
C Each plate is described by a set of 32 planar triangular panels
C  of equal size. The shell solution points are the centres of the 
C  panels. The solution is sought at the exterior points 
C  (0.05,0.05,z), the axis through the plates with z=0.001,0.003,
C  0.005, 0.007, 0.009.

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 LEBEM3: Subroutine for solving the exterior Laplace
C  equation (file LEBEM3.FOR contains the subroutine LEBEM3)
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  LSEM3T
      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=10)
C   Limit on the number of vertices (equal to the number of panels)
      INTEGER    MAXNV
      PARAMETER (MAXNV=12)
C   Limit on the number of points exterior to the plate, where 
C    potential properties are sought
      INTEGER    MAXNPE
      PARAMETER (MAXNPE=5)


C  Geometrical description of the plate(s)
C   Number of panels and counter
      INTEGER    NH,IPI
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    HELV(MAXNH,3)
C   The points exterior to the plate(s) where the potential 
C    properties 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
C    Coordinates of the exterior points
      REAL*8     PEXT(MAXNPE,3)
C  Shell Conditions
C   a(p)+b(p)=f(p)
C    function a at the central points
      REAL*8     HA(MAXNH)
C    function b at the central points
      REAL*8     HB(MAXNH)
C    function f at the central points
      REAL*8     HF(MAXNH)
C   A(p)+B(p)=F(p)
C    function A at the central points
      REAL*8     HAA(MAXNH)
C    function B at the central points
      REAL*8     HBB(MAXNH)
C    function F at the central points
      REAL*8     HFF(MAXNH)

C  Incident field
C    The incident potential {\phi} at the central points
      REAL*8     HIPHI(MAXNH)
C    The derivative of the incident potential {\phi} at the 
C     central points
      REAL*8     HIVEL(MAXNH)
C    The incident potential {\phi} at the exterior points
      REAL*8     PEIPHI(MAXNPE)

C Validation and control parameters
      LOGICAL    LSOL
      LOGICAL    LVALID
      REAL*8     EGEOM

C Solution
      REAL*8     PHIDIF(MAXNH)
      REAL*8     PHIAV(MAXNH)
      REAL*8     VELDIF(MAXNH)
      REAL*8     VELAV(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,z coordinates
      INTEGER    ICOORD
C  Counter through the exterior points
      INTEGER    IPE

C  Distance between plates
      REAL*8 GAP



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

C Describe the nodes and panels that make up the plate
C  :The two 1x1 plates are each divided into 32 panels, 
C  : (NH=64). VERTEX and HELV are defined anti-clockwise around 
C  : the plate when viewed from above so that the normal to the 
C  : plate is assumed to be upward
C  :Set up nodes
C  : Set NH, the number of panels
      NH=8
C  : Set NV, the number of vertices (equal to the number of panels)
      NV=10
C  : Set coordinates of the nodes


C  : Set up VERTEX, the coordinates of the vertices of the panels
      DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,10)

C   Plate 0 (z=0)
     * / 0.000D0, 0.000D0, 0.000D0,
     *   1.000D0, 0.000D0, 0.000D0,
     *   1.000D0, 1.000D0, 0.000D0,
     *   0.000D0, 1.000D0, 0.000D0,
     *   0.500D0, 0.500D0, 0.000D0,

C   Plate 1 (z=1)
     *   0.000D0, 0.000D0, 0.100D0,
     *   1.000D0, 0.000D0, 0.100D0,
     *   1.000D0, 1.000D0, 0.100D0,
     *   0.000D0, 1.000D0, 0.100D0,
     *   0.500D0, 0.500D0, 0.100D0/


        
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  :  plate into the potential domain.
      DATA ((HELV(IPI,ICOORD),ICOORD=1,3),IPI=1,8)

C Plate 0 (z=0)
     * / 1, 2, 5,     2, 3, 5,    3, 4, 5,     4, 1, 5,
     
C Plate 1 (z=0.5)
     *   6, 7, 10,    7, 8, 10,   8, 9,10,     9, 6, 10/

C Set NPE
      NPE=5
C Set coordinates of chosen exterior points
      GAP=0.1D0
      DO 300 IPE=1,NPE
        PEXT(IPE,1)=0.5D0
        PEXT(IPE,2)=0.5D0
        PEXT(IPE,3)=GAP*DFLOAT(2*IPE-1)/DFLOAT(2*NPE)
300   CONTINUE     

     
C Set the potential on the plates
      DO 310 IPI=1,NH
        HA(IPI)=1.0D0
        HB(IPI)=0.0D0
        HF(IPI)=0.0D0
        HAA(IPI)=1.0D0
        HBB(IPI)=0.0D0
C  Plate 1
        IF (IPI.LE.4) HFF(IPI)=0.0D0
C  Plate 2
        IF (IPI.GT.4) HFF(IPI)=1.0D0
310   CONTINUE

      DO 320 IPI=1,NH
        HIPHI(IPI)=0.0D0
        HIVEL(IPI)=0.0D0
320   CONTINUE     
      DO 330 IPE=1,NPE
        PEIPHI(IPE)=0.0D0
330   CONTINUE


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

      CALL LSEM3(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)

999   FORMAT(3F14.4)
      OPEN(10,FILE='LSEM3.OUT')
      WRITE(10,*) 'LSEM3: Test Problem Results'
      WRITE(10,*) '         z        computed pot    exact/app pot'
      DO 600 IPI=1,NPE
        WRITE(10,999) PEXT(IPI,3),PEPHI(IPI),PEXT(IPI,3)*10
600   CONTINUE



      END



                                                   