C*************************************************************** C Subroutine LSEMA 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/LSEMA.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 subroutine computes the solution of Laplace's equation C __ 2 C \/ {\phi} = 0 C C in the region exterior to a set of thin surfaces. C C C How to use the subroutine C ------------------------- C C The following diagram shows how the subroutine is to be used. C C .................................... C : : C : : C ---------------------- : -------------------------- : C | | : | | : C | MAIN PROGRAM |------->-----| LSEMA | : C | (e.g. LSEMA_t.for) | : | | : C | | : -------------------------- : C ---------------------- : | : C : > : C : | : C : ------------------------ : C Package ---------------->: | subordinate routines | : C : ------------------------ : C : : C : (this file) : C :..................................: C / | | C |_ > > C / | | C ................ ................ ................ C : : : -------- : : --------- : C : (geom3d.for) :---<---: |L3ALC | : : | GLS | : C : (geom2d.for) : : -------- : : --------- : C :..............: : -------------: : -------------: C : |subordinate|: : |subordinate|: C : | routines |: : | routines |: C : -------------: : -------------: C : : : : C : (l3alc.for) : : (gls.for) : C :..............: :..............: C C C The contents of the main program must be linked to LSEMA.FOR, L3ALC.FOR C GLS2.FOR, GL8.FOR, GEOM3D.FOR and GEOM2D.FOR. C C Method of solution C ------------------ C C In the main program, the shell(s) must be described as a set of C panels. The panels are defined by two indices (integers) which C label a node or vertex on the shell generator. The data structure C VERTEX lists and enumerates the coordinates of the vertices, the data C structure HELV defines each panel by indicating the labels for C the three nodes that are its vertices and hence enumerates the C panels. C C Format and parameters of the subroutine C --------------------------------------- C C The subroutine has the form C C SUBROUTINE LSEMA(MAXNV,NV,VERTEX,MAXNH,NH,HELV, C * MAXNPE,NPE,PEXT, C * HA,HB,HF,HAA,HBB,HFF, C * HIPHI,HIVEL,PEIPHI, C * LSOL,LVALID,EGEOM, C * PHIDIF,PHIAV,VELDIF,VELAV,PEPHI, C * AMAT,BMAT,L_EH, M_EH, C * PERM, XORY,C,WKSPC1,WKSPC2,WKSPC3) C The parameters to the subroutine C ================================ C Geometry of the shell {\Pi} (input) C integer MAXNV: The limit on the number of vertices of the triangles C that defines (approximates) {\Pi}. MAXNV>=2. C integer NV: The number of vertices on {\Pi}. 2<=NV<=MAXNV. C real VERTEX(MAXNV,2): The coordinates of the vertices. VERTEX(i,1), C VERTEX(i,2) are the x,y,z coordinates of the i-th vertex. C integer MAXNH: The limit on the number of panels describing {\Pi}. C MAXNH>=1. C integer NH: The number of panels describing {\Pi}. 1<=NH<=MAXNH. C integer HELV(MAXNH,2): The indices of the two vertices on the C generator defining each panel. The i-th panel has vertices C (VERTEX(HELV(i,1),1),VERTEX(HELV(i,1),2))) and C (VERTEX(HELV(i,2),1),VERTEX(HELV(i,2),2))). C Exterior points at which the solution is to be observed (input) C integer MAXNPE: Limit on the number of points exterior to the C shell. MAXNPE>=1. C integer NPE: The number of exterior points. 0<=NPE<=MAXNPE. C real PEXT(MAXNPE,2). The coordinates of the exterior point. C PEXT(i,1),PEXT(i,2) are the x,y,z coordinates of the i-th point. C Shell Conditions C a(p){\delta}(p)+b(p){\nu}(p)=f(p) C real HA(MAXNH). Function a at the central points C real HB(MAXNH). Function b at the central points C real HF(MAXNH). Function f at the central points C A(p){\Phi}(p)+B(p)V(p)=F(p) C real HAA(MAXNH). Function A at the central points C real HBB(MAXNH). Function B at the central points C real HFF(MAXNH). Function F at the central points C Incident field C real HIPHI(MAXNH). The incident potential {\phi} at the central C points C real HIVEL(MAXNH). The derivative of the incident potential {\phi} C at the central points C real PEIPHI(MAXNPE). The incident potential {\phi} at the exterior C points C Validation and control parameters (input) C logical LSOL: A switch to control whether the particular solution is C required. C logical LVALID: A switch to enable the choice of checking of C subroutine parameters. C real EGEOM: The maximum absolute error in the parameters that C describe the geometry. C Solution (output) C real SPHI(MAXNH): The potential ({\phi}) at the centres of the C boundary panels. C real SVEL(MAXNH): The (v or d{\phi}/dn where n is the outward C normal to the boundary) at the centres of the boundary C panels. C real PEPHI(MAXNPE): The potential ({\phi}) at the exterior points. C Working space C REAL*8 AMAT(2*MAXNH,2*MAXNH) C REAL*8 BMAT(2*MAXNH,2*MAXNH) C REAL*8 L_EH(MAXNPE,MAXNH) C REAL*8 M_EH(MAXNPE,MAXNH) C REAL*8 C(2*MAXNH) C REAL*8 HALPHA(2*MAXNH) C REAL*8 HBETA(2*MAXNH) C REAL*8 HFFF(2*MAXNH) C REAL*8 PHIINF(2*MAXNH) C REAL*8 VELINF(2*MAXNH) C REAL*8 WKSPC1(2*MAXNH) C REAL*8 WKSPC2(2*MAXNH) C REAL*8 WKSPC3(2*MAXNH) C Notes on the geometric parameters C --------------------------------- C (1) Each of the vertices listed in VERTEX must be distinct points C with respect to EGEOM. C (2) The boundary must be complete and closed. Thus C HELV(i,2)=HELV(i+1,1) for i=1..NH-1 and VERTEX(HELV(1,1),1)=0 C and VERTEX(HELV(NH,2),1)=0. C (3) The indices of the nodes listed in HELV must be such that they C are ordered counter-clockwise around the generator of the boundary. C (4) The generator of the largest panel must be no more than 10x C the length of the generator of the smallest panel. C Notes on the exterior points C ---------------------------- C (1) The points in PEXT should lie outside the shell, as defined C by the parameters VERTEX and HELV. Any point lying outside the C shell. C The subroutine SUBROUTINE LSEMA(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) PARAMETER (MAXNGQ=32) PARAMETER (MAXNTQ=10000) C Boundary geometry C Limit on the number of vertices on {\Pi} INTEGER MAXNV C The number of vertices on {\Pi} INTEGER NV C The coordinates of the vertices on {\Pi} REAL*8 VERTEX(MAXNV,2) C Limit on the number of panels describing {\Pi} INTEGER MAXNH C The number of panels describing {\Pi} INTEGER NH C The indices of the vertices describing each panel INTEGER HELV(MAXNH,2) C Exterior points at which the solution is to be observed C Limit on the number of points exterior to the shell where C solution is sought INTEGER MAXNPE C The number of exterior points INTEGER NPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,2) C Shell Conditions C a(p){\delta}(p)+b(p){\nu}(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){\Phi}(p)+B(p)V(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) REAL*8 C(2*MAXNH) REAL*8 HALPHA(2*MAXNH) REAL*8 HBETA(2*MAXNH) REAL*8 HFFF(2*MAXNH) REAL*8 PHIINF(2*MAXNH) REAL*8 VELINF(2*MAXNH) REAL*8 WKSPC1(2*MAXNH) REAL*8 WKSPC2(2*MAXNH) REAL*8 WKSPC3(2*MAXNH) C Further information from system solution GLS INTEGER*4 PERM(2*MAXNH) LOGICAL XORY(2*MAXNH) REAL*8 DIST2 REAL*8 DIAM C Constants C Real scalars: 0, 1, 2, half, pi REAL*8 ZERO,ONE,TWO,PI C Geometrical description of the shell C panels counter INTEGER IPI,JPI C The points exterior to the shell where the solution is sought INTEGER IPE C Parameters for L3ALC REAL*8 P(2),PA(2),PB(2),QA(2),QB(2),QC(2),VECP(2) LOGICAL LPONEL C Quadrature rule information C [Note that in this program two quadrature rules are used: one for C the case when the point P lies on the panel (LPONEL=.TRUE.) and C one for the case when P does not lie on the panel.] C Quadrature rule used when LPONEL=.TRUE. C Number of quadrature points INTEGER NGQON C Abscissae of the actual quadrature rule REAL*8 AGQON(MAXNGQ) C Weights of the actual quadrature rule REAL*8 WGQON(MAXNGQ) C Quadrature rule used when LPONEL=.FALSE. C Number of quadrature points INTEGER NGQOFF C Abscissae of the actual quadrature rule REAL*8 AGQOFF(MAXNGQ) C Weights of the actual quadrature rule REAL*8 WGQOFF(MAXNGQ) C Quadrature rule parameters for L3ALC C Actual number of quadrature points INTEGER NGQ C Abscissae of the actual quadrature rule in the generator direction REAL*8 AGQ(MAXNGQ) C Weights of the actual quadrature rule REAL*8 WGQ(MAXNGQ) C Counter through the quadrature points INTEGER IGQ C Abscissae of the actual quadrature rule in the theta direction REAL*8 ATQ(MAXNTQ) C Weights of the actual quadrature rule REAL*8 WTQ(MAXNTQ) C Validation and control parameters for subroutine L3ALC LOGICAL LVAL REAL*8 EQRULE LOGICAL LL LOGICAL LM LOGICAL LMT LOGICAL LN C Parameters for subroutine L3ALC. REAL*8 DISL REAL*8 DISM REAL*8 DISMT REAL*8 DISN C Other variables C Failure flag LOGICAL LFAIL C Accumulation of solution {\phi} REAL*8 SUMPHI C Error flag LOGICAL LERROR REAL*8 WKSPCE(2*MAXNTQ+MAXNGQ) REAL*8 SIZMAX,SIZMIN C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 PI=4.0D0*ATAN(ONE) C Set up validation and control parameters C Switch off the validation of L3ALC LVAL=.FALSE. C Set EQRULE EQRULE=1.0D-6 C Set EGEOM EGEOM=1.0D-6 C Set up the quadrature rule(s). C Set up quadrature rule for the case when P is not on the panel. C Set up 8 point Gauss-Legendre rules CALL GL8(MAXNGQ,NGQOFF,WGQOFF,AGQOFF) C Set up quadrature rule for the case when P is on the panel. C This is done by splitting the quadrature rule at the centre. NGQON=2*NGQOFF DO 330 IGQ=1,NGQOFF AGQON(IGQ)=AGQOFF(IGQ)/TWO AGQON(NGQOFF+IGQ)=0.5D0+AGQOFF(IGQ)/TWO WGQON(IGQ)=WGQOFF(IGQ)/TWO WGQON(NGQOFF+IGQ)=WGQOFF(IGQ)/TWO 330 CONTINUE C Validation C ========== C Validation of parameters of LEBEM3 C --------------------------------- IF (LVALID) THEN C Validation of main paramters LERROR=.FALSE. IF (MAXNV.LT.2) THEN WRITE(*,*) 'MAXNV = ',MAXNV WRITE(*,*) 'ERROR(LSEMA) - must have MAXNV>=4' LERROR=.TRUE. END IF IF (NV.LT.2.OR.NV.GT.MAXNV) THEN WRITE(*,*) 'NV = ',NV WRITE(*,*) 'ERROR(LSEMA) - must have 4<=NV<=MAXNV' LERROR=.TRUE. END IF IF (MAXNH.LT.2) THEN WRITE(*,*) 'MAXNH = ',MAXNH WRITE(*,*) 'ERROR(LSEMA) - must have MAXNH>=2' LERROR=.TRUE. END IF IF (NH.LT.3.OR.NH.GT.MAXNH) THEN WRITE(*,*) 'NH = ',NH WRITE(*,*) 'ERROR(LSEMA) - must have 4<=NH<=MAXNH' LERROR=.TRUE. END IF IF (MAXNPE.LT.1) THEN WRITE(*,*) 'MAXNPE = ',MAXNPE WRITE(*,*) 'ERROR(LSEMA) - must have MAXNPE>=1' LERROR=.TRUE. END IF IF (NPE.LT.0.OR.NPE.GT.MAXNPE) THEN WRITE(*,*) 'NPE = ',NPE WRITE(*,*) 'ERROR(LSEMA) - must have 3<=NPE<=MAXNPE' LERROR=.TRUE. END IF IF (EGEOM.LE.ZERO) THEN WRITE(*,*) 'NPE = ',NPE WRITE(*,*) 'ERROR(LSEMA) - EGEOM must be positive' LERROR=.TRUE. END IF IF (LERROR) THEN LFAIL=.TRUE. WRITE(*,*) WRITE(*,*) 'Error(s) found in the main parameters of LSEMA' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Find the diameter DIAM of the boundary DIAM=0.0 DO 100 IV=1,NV-1 PA(1)=VERTEX(IV,1) PA(2)=VERTEX(IV,2) DO 110 JV=IV+1,NV PB(1)=VERTEX(JV,1) PB(2)=VERTEX(JV,2) DIAM=MAX(DIAM,DIST2(PA,PB)) 110 CONTINUE 100 CONTINUE IF (LVALID) THEN LERROR=.FALSE. C Check that the boundary defined by HELV is complete and closed DO 120 IPI=1,NH 120 CONTINUE C Check that EGEOM is not too large IF (EGEOM.GT.DIAM/100.0D0) THEN WRITE(*,*) 'EGEOM = ',EGEOM WRITE(*,*) 'ERROR(LSEMA) - EGEOM is set too large' LERROR=.TRUE. END IF IF (LERROR) THEN LFAIL=.TRUE. WRITE(*,*) WRITE(*,*) 'Error in boundary geometry or EGEOM' WRITE(*,*) 'Execution terminated' END IF END IF IF (LVALID) THEN C Check that the vertices are distinct with respect to EGEOM LERROR=.FALSE. DO 130 IV=1,NV-1 PA(1)=VERTEX(IV,1) PA(2)=VERTEX(IV,2) DO 140 JV=IV+1,NV PB(1)=VERTEX(JV,1) PB(2)=VERTEX(JV,2) IF (ABS(PA(1)-PB(1)).LT.EGEOM) THEN IF (ABS(PA(2)-PB(2)).LT.EGEOM) THEN WRITE(*,*) 'Vertices ',IV,JV,' are not distinct' LERROR=.TRUE. END IF END IF 140 CONTINUE 130 CONTINUE IF (LERROR) THEN WRITE(*,*) WRITE(*,*) 'ERROR(LSEMA) - Vertices (see above) coincide' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Check that the panels are not of disproportionate sizes IF (LVALID) THEN SIZMAX=ZERO SIZMIN=DIAM DO 150 IPI=1,NH QA(1)=VERTEX(HELV(IPI,1),1) QA(2)=VERTEX(HELV(IPI,1),2) QB(1)=VERTEX(HELV(IPI,2),1) QB(2)=VERTEX(HELV(IPI,2),2) SIZMAX=MAX(SIZMAX,DIST2(QA,QB)) SIZMIN=MIN(SIZMIN,DIST2(QA,QB)) 150 CONTINUE IF (SIZMAX.GT.10.0D0*SIZMIN) THEN WRITE(*,*) 'WARNING(LSEMA) - panels of disproportionate' WRITE(*,*) ' sizes' END IF END IF C Check that the boundary does not contain sharp angles IF (LVALID) THEN LERROR=.FALSE. DO 160 IPI=2,NH IF (HELV(IPI-1,2).EQ.HELV(IPI,1)) THEN QA(1)=VERTEX(HELV(IPI,1),1) QA(2)=VERTEX(HELV(IPI,1),2) QB(1)=VERTEX(HELV(IPI,2),1) QB(2)=VERTEX(HELV(IPI,2),2) QC(1)=VERTEX(HELV(IPI-1,1),1) QC(2)=VERTEX(HELV(IPI-1,1),2) IF (DIST2(QC,QB).LT.MAX(DIST2(QA,QB),DIST2(QA,QC))) THEN WRITE(*,*) 'Sharp angle at node ',HELV(IPI,1) LERROR=.TRUE. END IF END IF 160 CONTINUE IF (LERROR) THEN WRITE(*,*) WRITE(*,*) 'WARNING(LSEMA) - Boundary has sharp angles' END IF END IF C Validation of the surface functions IF (LVALID.AND.LSOL) THEN LERROR=.FALSE. DO 170 IPI=1,NH IF (MAX(ABS(HA(IPI)),ABS(HB(IPI))).LT.1.0D-6) * LERROR=.TRUE. 170 CONTINUE IF (LERROR) THEN WRITE(*,*) WRITE(*,*) 'ERROR(LSEMA) - at most one of HA(i),HB(i)' WRITE(*,*) ' may be zero for all i' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Compute the discrete L, M, Mt and N matrices C Loop(IPI) through the points on the boundary DO 510 IPI=1,NH C Set P PA(1)=VERTEX(HELV(IPI,1),1) PA(2)=VERTEX(HELV(IPI,1),2) PB(1)=VERTEX(HELV(IPI,2),1) PB(2)=VERTEX(HELV(IPI,2),2) P(1)=(PA(1)+PB(1))/TWO P(2)=(PA(2)+PB(2))/TWO C Set VECP to the normal on the boundary of the panel at P CALL NORM2(PA,PB,VECP) C Loop(IPI) through the panels DO 520 JPI=1,NH C Set QA and QB, the coordinates of the edges of the JSEth panel QA(1)=VERTEX(HELV(JPI,1),1) QA(2)=VERTEX(HELV(JPI,1),2) QB(1)=VERTEX(HELV(JPI,2),1) QB(2)=VERTEX(HELV(JPI,2),2) C Set LPONEL IF (IPI.EQ.JPI) THEN LPONEL=.TRUE. ELSE LPONEL=.FALSE. END IF C Select quadrature rule for L3ALC C : Select the quadrature rule AGQON-WGQON in the case when the C : point p lies on the panel, otherwise select AGQOFF-WGQOFF C [Note that the overall method would benefit from selecting from C a wider set of quadrature rules, and an appropriate method C of selection] IF (LPONEL) THEN NGQ=NGQON DO 600 IGQ=1,NGQ AGQ(IGQ)=AGQON(IGQ) WGQ(IGQ)=WGQON(IGQ) 600 CONTINUE ELSE NGQ=NGQOFF DO 610 IGQ=1,NGQ AGQ(IGQ)=AGQOFF(IGQ) WGQ(IGQ)=WGQOFF(IGQ) 610 CONTINUE END IF C Quadrature rule in the theta direction is constructed out of individual C Gauss rules so that the length of each is approximately equal to the C length of the panel at the generator. RADMID=(QA(1)+QB(1))/TWO SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+ * (QA(2)-QB(2))*(QA(2)-QB(2)) GLEN=SQRT(SGLEN) CIRMID=PI*RADMID NDIV=1+CIRMID/GLEN TDIV=ONE/DBLE(NDIV) NTQ=NDIV*NGQ IF (NTQ.GT.MAXNTQ) THEN WRITE(*,*) 'ERROR(LSEMA) - MAXNTQ is set too small' STOP END IF DO 145 IDIV=1,NDIV DO 155 IGQ=1,NGQ WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV) ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+ * TDIV*DBLE(IDIV-1) 155 CONTINUE 145 CONTINUE C All operators are required LL=.TRUE. LM=.TRUE. LMT=.TRUE. LN=.TRUE. C Call of L3ALC routine to compute [L], [M], [Mt], [N] CALL L3ALC(P,VECP,QA,QB,LPONEL, * MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ, * LVAL,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN, * WKSPCE) AMAT(IPI,JPI)=-DISM AMAT(IPI,NH+JPI)=0.0D0 AMAT(NH+IPI,JPI)=-DISN AMAT(NH+IPI,NH+JPI)=0.0D0 BMAT(IPI,JPI)=-DISL BMAT(IPI,NH+JPI)=0.0D0 BMAT(NH+IPI,JPI)=-DISMT BMAT(NH+IPI,NH+JPI)=0.0D0 C Close loop(JPI) 520 CONTINUE AMAT(IPI,NH+IPI)=1.0D0 BMAT(NH+IPI,NH+IPI)=-1.0D0 C Close loop(IPI) 510 CONTINUE C Set incident field to zero DO 550 IPI=1,NH HIPHI(IPI)=0.0D0 HIVEL(IPI)=0.0D0 550 CONTINUE DO 560 IPI=1,NPE PEIPHI(IPI)=0.0D0 560 CONTINUE C Set C,HALPHA,HBETA,HFFF DO 700 IPI=1,NH C Set C as the incident potential C(IPI)=HIPHI(IPI) C(NH+IPI)=HIVEL(IPI) C Set HALPHA as the coefficient of [{\delta} {\Phi}] HALPHA(IPI)=HA(IPI) HALPHA(NH+IPI)=HAA(IPI) C Set HBETA as the coefficient of [{\nu} V] HBETA(IPI)=HB(IPI) HBETA(NH+IPI)=HBB(IPI) C Set HFFF as [f F] HFFF(IPI)=HF(IPI) HFFF(NH+IPI)=HFF(IPI) 700 CONTINUE C On solution C becomes {\zeta}/{\epsilon} at the central points CALL GLS(2*MAXNH,2*NH,AMAT,BMAT,C,HALPHA,HBETA, * HFFF,PHIINF,VELINF,LFAIL,PERM,XORY,WKSPC1,WKSPC2,WKSPC3) C Set C DO 710 IPI=1,NH PHIDIF(IPI)=PHIINF(IPI) PHIAV(IPI)=PHIINF(NH+IPI) VELDIF(IPI)=VELINF(IPI) VELAV(IPI)=VELINF(NH+IPI) 710 CONTINUE C Compute sound pressures at the selected exterior points. C Loop through the the points in the exterior region DO 800 IPE=1,NPE C Set P P(1)=PEXT(IPE,1) P(2)=PEXT(IPE,2) C Set VECP, this is arbitrary as the derivative of the potential at P C is not sought. VECP(1)=ONE VECP(2)=ZERO C Initialise SUMPHI to the incident potential SUMPHI=PEIPHI(IPE) C Loop(JPI) through the panels DO 850 JPI=1,NH C Compute the discrete L and M integral operators. C Set QA and QB, the coordinates of the edges of the JSEth panel QA(1)=VERTEX(HELV(JPI,1),1) QA(2)=VERTEX(HELV(JPI,1),2) QB(1)=VERTEX(HELV(JPI,2),1) QB(2)=VERTEX(HELV(JPI,2),2) C Quadrature rule in the generator direction NGQ=NGQOFF DO 1610 IGQ=1,NGQ AGQ(IGQ)=AGQOFF(IGQ) WGQ(IGQ)=WGQOFF(IGQ) 1610 CONTINUE C Quadrature rule in the theta direction is constructed out of individual C Gauss rules so that the length of each is approximately equal to the C length of the panel at the generator. RADMID=(QA(1)+QB(1))/TWO SGLEN=(QA(1)-QB(1))*(QA(1)-QB(1))+ * (QA(2)-QB(2))*(QA(2)-QB(2)) GLEN=SQRT(SGLEN) CIRMID=PI*RADMID NDIV=1+CIRMID/GLEN TDIV=ONE/DBLE(NDIV) NTQ=NDIV*NGQ IF (NTQ.GT.MAXNTQ) THEN WRITE(*,*) 'ERROR(LSEMA) - MAXNTQ is set too small' STOP END IF DO 1145 IDIV=1,NDIV DO 1155 IGQ=1,NGQ WTQ((IDIV-1)*NGQ+IGQ)=WGQ(IGQ)/DBLE(NDIV) ATQ((IDIV-1)*NGQ+IGQ)=AGQ(IGQ)/DBLE(NDIV)+ * TDIV*DBLE(IDIV-1) 1155 CONTINUE 1145 CONTINUE C All the points do not lie on the boundary hence LPONEL=.FALSE. LPONEL=.FALSE. C Only L, M operators are required. Set LL,LM true, C LMT,LN false. LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Call L3ALC. CALL L3ALC(P,VECP,QA,QB,LPONEL, * MAXNGQ,NGQ,AGQ,WGQ,MAXNTQ,NTQ,ATQ,WTQ, * LVAL,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN, * WKSPCE) C Accumulate phi SUMPHI=SUMPHI+DISM*PHIDIF(JPI)-DISL*VELDIF(JPI) L_EH(IPE,JPI)=DISL M_EH(IPE,JPI)=DISM C Close loop (JPI) through the panels 850 CONTINUE PEPHI(IPE)=SUMPHI C Close loop(IPE) through the exterior points 800 CONTINUE END C ---------------------------------------------------------------------- C Subordinate routines for LEBEMA C =============================== C Subroutines required for L3ALC (not in file L3ALC.FOR) C Subroutine for returning the square root. REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END