program main
c*********************************************************************72
c
cc MAIN is the main program for the BETIS boundary element code.
c
c Discussion:
c
c BETIS is a program for the solution of Laplace's equation,
c using the boundary element method.
c
c BETIS was developed by Federico Paris and Jose Canas, of
c the department of elasticity and strength of materials,
c in the Industrial Engineering school of the University of
c Seville, Spain.
c
c The program uses linear continuous elements and Dirichlet, Neumann
c as well as mixed boundary conditions can be taken into
c consideration.
c
c The original version of BETIS was written in 1978.
c
c Copies of BETIS may be found at
c
c http://www.esi2.us.es/mmc/erm
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7,
c LC: TA347.B69.P34.
c
c Parameters:
c
c Local, real A(2), stores the results of the
c integrations of the kernels along an element.
c
c Local, real B(2), stores the results of the
c integrations of the kernels along an element.
c
c Local, real CARGA(NPMAX), the load vector of the system of equations
c
c Local, real COEF(NPMAX,NPMAX), the coefficients of the final system of
c equations to be solved
c
c Local, real DFI(NPMAX,2), the fluxes at the nodes of the boundary. It has
c double the dimension of FI because two fluxes at each node
c are stored.
c
c Local, real FI(NPMAX), the potential values at the nodes of the boundary.
c
c Local, real FINT(NPMAX-1), values of the potential at the internal points.
c
c User input, character*60 IENC, the title of the problem.
c
c User input, integer IGAUSS, the number of points to be used
c in the Gauss quadrature rule. A value of 4 is typical.
c 1 <= IGAUSS <= 50 is required.
c
c Local, integer IMP, the unit number for the output data file.
c
c Local, integer IPIVO(NPMAX), workspace used for pivoting during
c the solution of the linear system.
c
c Local, integer LEC, the unit number for the input data file.
c
c User input, integer NANU, the index of the node where the potential is
c going to be considered zero. This is only for the Neumann problem, the
c fluxes along the boundary being specified.
c
c Local, integer NCOD(NPMAX), the codes of the nodes of the boundary used to
c identify the boundary conditions.
c
c User input, integer NP, the number of points (nodes) of the boundary where
c the boundary integral equation is going to be applied. This number
c represents also the number of elements used in the discretization of the
c boundary.
c
c User input, integer NPI, the number of internal points where the
c potential is going to be calculated.
c
c Local, integer parameter NPMAX, the maximum number of elements (plus one)
c that can be used to model the boundary of a certain problem.
c
c Local, real X(*), the X coordinates of the boundary points
c used in the discretization.
c
c Local, real XINT(NPMAX-1), the X coordinates of the internal points where
c the potential is going to be calculated
c
c Local, real Y(NPMAX), the Y coordinates of the boundary points
c used in the discretization.
c
c Local, real YINT(NPMAX-1), the Y coordinates of the internal points where
c the potential is going to be calculated
c
implicit none
integer npmax
parameter ( npmax = 101 )
real a(2)
real ax
real ay
real b(2)
real carga(npmax)
real cdiag
real coef(npmax,npmax)
real cose
real den
real dfi(npmax,2)
real dfx
real dfy
real fi(npmax)
real fint(npmax-1)
integer i
integer ielem
character*60 ienc
integer igauss
integer imp
integer ipivo(npmax)
integer j
integer jj
integer k
integer kk
integer kn
integer knma
integer knme
integer l
integer lec
integer na
character*15 namedat
character*15 namesal
integer nanu
integer nc
integer ncod(npmax)
integer nint
integer nodo
integer np
integer npi
real omeg(50)
real x(npmax)
real xi(50)
real xint(npmax-1)
real xnue
real y(npmax)
real yint(npmax-1)
real ynue
lec=10
imp=11
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BETIS'
write ( *, '(a)' ) ' FORTRAN77 version'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' A 2D boundary element program'
write ( *, '(a)' ) ' by Federico Paris and Jose Canas.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'Data File name?'
read(*,'(a)') namedat
open(lec,file=namedat,status='old',err=9988)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'Result File name?'
read(*,'(a)') namesal
open(imp,file=namesal,status='new')
write ( *, '(a)' ) 'Number of Gauss Points?'
read(*,*) igauss
call gauss_qn ( igauss, xi, omeg )
c
c Reading of the basic variables of the program
c
c The title of the example is written (maximum 50 characters) in the
c first line of the data file and read as the variable ienc
c
read(lec,'(a)') ienc
c
c The second line of the data file includes the number of nodes of
c the boundary (np), the number of internal points where the
c potential is going to be calculated (npi) and the node where the
c potential is going to be forced to be zero (nanu).
c
c np cannot be zero. npi is zero when no information about
c the domain is required. nanu is zero when the position of the
c potential distribution is not undetermined, which is the case
c of mixed or Dirichlet boundary conditions. In the case of Neumann
c boundary conditions the variable must include the number of the
c node where the potential is going to be taken as null.
c
c The three variables must be written as integers separated by a
c blank space
c
read(lec,*) np,npi,nanu
write(imp,1010) ienc
write(imp,1001)np,npi,igauss,nanu
c
c The main arrays are initialized to zero
c
do i = 1, np
fi(i) = 0.0E+00
dfi(i,1) = 0.0E+00
dfi(i,2) = 0.0E+00
end do
c
c The coordinates of the internal points are read, if it has
c previously been specified that this is required (npi not equal
c zero).
c
c As many lines as indicated by the variable npi must necessarily
c be included. Each of these lines must have the coordinates x and y
c of each internal point in real format separated by a blank space.
c
if ( npi .ne. 0 ) then
read(lec,*) (xint(i),yint(i),i=1,npi)
end if
c
c Reading and generation of the coordinates of the nodes of the
c boundary.
c
c Each line of this block must have the number of the node (integer)
c and its two coordinates x and y (real), separated by a blank space.
c
c The nodes of the boundary must be numbered anticlockwise,
c starting at any point. The numbers of the nodes in the subsequent
c lines must always increase although they are not required to be
c consecutive.
c If two lines have non consecutive node numbers the program
c will automatically generate the coordinates of the intermediate
c points dividing the straight line between the two specified nodes
c into segments of equal length. The first node can also be used for
c automatic generation of the last nodes of the boundary, giving it
c the number np+1 and being obviously the last line of this block.
c
l=0
2 continue
read(lec,*)i,x(i),y(i)
if ( 0 < l ) then
nint = i - l
ax = ( x(i) - x(l) ) / nint
ay = ( y(i) - y(l) ) / nint
end if
4 continue
l = l + 1
if(i-l) 5,6,7
7 x(l) = x(l-1)+ax
y(l) = y(l-1)+ay
go to 4
6 if(np-i) 8,9,2
9 x(np+1) = x(1)
y(np+1) = y(1)
8 continue
c
c Reading and generation of the boundary conditions of the problem
c
c Each line of this block must include the number of the node
c (integer), the code of the node (integer), the potential (real),
c the flux before(real) and the flux after the node (real). Never, of
c course, will the three real variables be known. When a variable is
c unknown a zero value is placed, the code identifying which is/are
c the known variable/s.
c The following list will help with the preparation of the data.
c
c Code Variable/s to be specified
c
c 1 flux before the node [dfi(i,1)], flux after [dfi(i,2)]
c
c 2 potential [fi(i)], flux after the node[dfi(i,2)]
c
c 3 potential [fi(i)], flux before the node[dfi(i,1)]
c
c 4 potential [fi(i)], the external normal is discontinuous at i
c
c 5 potential [fi(i)], the external normal is continuous at i
c
c A linear automatic generation similar to the one used for the nodes
c coordinates is also employed here for the boundary conditions. It
c can be used independently of the geometry, because the values of
c the variables at the intermediate nodes are generated using the
c number of the nodes and not with its geometrical position.
c
c The same general features explained for the generation of the
c coordinates are applicable in this case.
c
l = 0
10 continue
read(lec,*) i,ncod(i),fi(i),dfi(i,1),dfi(i,2)
if(i-l-1)11,12,13
13 nint = i - l
ax = ( fi(i) - fi(l) ) / nint
ay = ( dfi(i,1) - dfi(l,2)) / nint
12 l = l + 1
if(i-l)11,14,15
15 if ( ncod(i) .le. 3 .and. ncod(i) .ne. 2 ) go to 16
ncod(l) = 5
fi(l) = fi(l-1) + ax
go to 12
16 ncod(l) = 1
dfi(l,1) = dfi(l-1,2) + ay
dfi(l,2) = dfi(l,1)
go to 12
14 continue
if(np-i)17,18,10
c
c The first node is also used as the last plus one to simplify later
c instructions.
c
18 fi(np+1) = fi(1)
dfi(np+1,1) = dfi(1,1)
dfi(np+1,2) = dfi(1,2)
ncod(np+1) = ncod(1)
17 continue
c
c Writing of input data in order to check them.
c
19 write(imp,1002)
write(imp,1003) (i,x(i),y(i),i=1,np)
write(imp,1102)
write(imp,1103) (i,ncod(i),fi(i),(dfi(i,k),k=1,2),i=1,np)
c
c Starting of the application of the Boundary Element Method.
c
c The variable NODO represents the node of the boundary where the
c integral equation is being generated.
c
do nodo = 1, np
c
c The row of the coefficient matrix corresponding to the equation
c to be generated and the corresponding position in the load vector
c are initialized.
c
c The variable CDIAG, also initialized, is used to store the sum of
c the A coefficients corresponding to the point where the integral
c equation is being applied. In this way, the calculation of the free
c term is avoided.
c
do i = 1, np
coef(nodo,i) = 0.0E+00
end do
cdiag = 0.0E+00
carga(nodo) = 0.0E+00
c
c The variable IELEM represents the element along which the
c integration is going to be performed, from the node NODO.
c
do ielem = 1, np
if(nodo.eq.1.and.ielem.eq.np) go to 106
if(ielem.eq.nodo)go to 105
if((ielem+1).eq.nodo)go to 106
c
c The node where the integral equation is applied does not belong
c to the element along which the integration is performed, the
c integrations being performed numerically in the subroutine HGNUM.
c
call hgnum ( x(nodo), y(nodo), x(ielem), y(ielem),
& x(ielem+1), y(ielem+1), a(1), a(2), b(1), b(2), igauss,
& xi, omeg )
go to 160
c
c The node where the integral equation is applied belongs to the
c element along which the integration is performed, an analytical
c expression being used in the subroutine HGANA.
c
105 call hgana ( x(nodo), y(nodo), x(ielem+1), y(ielem+1),
& a(1), a(2), b(1), b(2) )
go to 160
106 call hgana ( x(ielem), y(ielem), x(nodo), y(nodo), a(1),
& a(2), b(1), b(2) )
ax = b(1)
b(1) = b(2)
b(2) = ax
c
c The values of the coefficients A are stored in CDIAG in order to
c calculate the free term, eqn 3.5.13, Section 3.5.2.3.
c
160 cdiag = cdiag+a(1)+a(2)
c
c The values of the integrations performed are going to be introduced
c into the system of equations according to the codes of the nodes
c of the elements along which the integrations have been performed.
c
c The variable K is used here to refer to the initial (K=1) and final
c node (K=2) of the element along which the integration has just been
c performed.
c
do k = 1, 2
kn = ielem + k - 1
if ( kn .eq. np + 1 ) then
kn = 1
end if
c
c The system of equations is built up according to the boundary
c conditions which are reflected in the codes.
c
c Code 1: the values of the derivative of the function before and
c after the node are known, case 1, Section 3.5.3.
c
if ( ncod(kn) .eq. 1 ) then
carga(nodo) = carga(nodo)+b(k)*dfi(kn,3-k)
coef(nodo,kn) = coef(nodo,kn)+a(k)
c
c Code 2: the value of the potential at the node and its derivative
c after the node are known, case 2a, Section 3.5.3.
c
else if ( ncod(kn) .eq. 2 ) then
if ( k .eq. 2 ) then
carga(nodo) = carga(nodo) - a(k) * fi(kn)
coef(nodo,kn) = coef(nodo,kn) - b(k)
else
carga(nodo) = carga(nodo)+b(k)*dfi(kn,3-k)-a(k)*fi(kn)
end if
c
c Code 3: the value of the potential at the node and its derivative
c before the node are known, case 2b, Section 3.5.3.
c
else if ( ncod(kn) .eq. 3 ) then
if ( k .eq. 2 ) then
carga(nodo) = carga(nodo)+b(k)*dfi(kn,3-k)-a(k)*fi(kn)
else
carga(nodo) = carga(nodo)-a(k)*fi(kn)
coef(nodo,kn) = coef(nodo,kn)-b(k)
end if
c
c Code 4: the value of the function at the node is known the normal
c at the node being discontinuous. The approach of the gradient is
c employed, case 3b, Section 3.5.3.
c
else if ( ncod(kn) .eq. 4 ) then
den = sqrt ( ( x(ielem) - x(ielem+1) )**2
& + ( y(ielem) - y(ielem+1) )**2 )
xnue = (y(ielem+1)-y(ielem))/den
ynue = (x(ielem)-x(ielem+1))/den
knma = kn+1
knme = kn-1
if(kn.eq.1)knme = np
if(kn.eq.np)knma = 1
den = (x(knma)-x(kn))*(y(knme)-y(kn))-(x(knme)-x(kn))*
& (y(knma)-y(kn))
c
c DFX and DFY represent the components of the gradient according
c to a linear assumption of the evolution of the function.
c
dfx = ((y(kn)-y(knma))*fi(knme)+(y(knma)-y(knme))*fi(kn)+
& (y(knme)-y(kn))*fi(knma))/den
dfy = ((x(knma)-x(kn))*fi(knme)+(x(knme)-x(knma))*fi(kn)+
& (x(kn)-x(knme))*fi(knma))/den
c
c If the potential is constant in the neighborhood of the node under
c study the direction of the gradient is non-determined. One is
c arbitrarily chosen.
c
if ( dfx .ne. 0.0E+00 .or. dfy .ne. 0.0E+00 ) then
cose = (dfx*xnue+dfy*ynue)/sqrt(dfx**2+dfy**2)
else
cose = 0.7071E+00
end if
c
c The variable cose represents the constants D that appear in eqn
c 3.5.19. It is associated to the cosine of the angle that forms the
c gradient whose components are DFX and DFY, with the flux at the
c node of the element under consideration.
c The variable cose is stored in DFI in order to avoid recalculating
c after the solution of the system of equations.
c
kk = 3-k
dfi(kn,kk) = cose
coef(nodo,kn) = coef(nodo,kn)-b(k)*cose
carga(nodo) = carga(nodo)-a(k)*fi(kn)
c
c Code 5: the value of the function at the node is known, the normal
c being continuous at the node, case 3a, Section 3.5.3.
c Code 5 is programmed using the same sentences as code 4, the
c artificial variable COSE receiving the value 1.
c
else if ( ncod(kn) .eq. 5 ) then
cose = 1.0E+00
kk = 3-k
dfi(kn,kk) = cose
coef(nodo,kn) = coef(nodo,kn)-b(k)*cose
carga(nodo) = carga(nodo)-a(k)*fi(kn)
end if
end do
c
c End of the integrations along all the elements from the node nodo.
c
end do
c
c The value of the free term is placed in its correct position
c according to the code of the node.
c
if ( ncod(nodo) .ne. 1 ) then
carga(nodo) = carga(nodo) + cdiag * fi(nodo)
else
coef(nodo,nodo) = coef(nodo,nodo)-cdiag
end if
c
c End of the generation of the integral equations.
c
end do
c
c If some potential is specified along the boundary the system is
c ready to be solved. If all the derivatives are known along the
c boundary the value of the potential must be forced to be zero at
c one point (node nanu), equation 3.4.16, Section 3.4.3.
c
na = 0
if ( nanu .ne. 0 ) then
na = 1
do i = 1, np
coef(np+1,i) = 0.0E+00
end do
coef(np+1,nanu) = 1.0E+00
carga(np+1) = 0.0E+00
end if
c
c The system of equations is now ready to be solved.
c
270 continue
call pivo ( coef, np, carga, ipivo, npmax, imp, na )
c
c The vector carga stores the solution of the system of equations,
c which must now be placed in its correct position.
c
do nodo = 1, np
nc = ncod(nodo)
c
c Code 1: CARGA stores the potential
c
if ( nc .eq. 1 ) then
fi(nodo) = carga(nodo)
c
c Code 2 and Code 3: CARGA stores the value of the derivative
c of the function before and after the node respectively.
c
else if ( nc .eq. 2 .or. nc .eq. 3 ) then
nc = nc - 1
dfi(nodo,nc) = carga(nodo)
c
c Code 4: CARGA stores the value of the gradient. The derivatives are
c obtained multiplying the gradient by the cosine of the angle that
c it forms with the corresponding normals, eqns 3.5.18, Section
c 3.5.3. The variable COSE was stored in the variable DFI.
c
c If the gradient forms 90 degrees with the normal the flux is zero
c and coincides with the value of the cosine previously stored.
c
else if ( nc .eq. 4 ) then
do i=1,2
if ( 1.0E-10 .lt. abs ( dfi(nodo,i) ) ) then
dfi(nodo,i) = carga(nodo) * dfi(nodo,i)
end if
end do
c
c Code 5: CARGA stores the value of the derivative of the potential
c before and after the node.
c
else if ( nc .eq. 5 ) then
do i = 1, 2
dfi(nodo,i) = carga(nodo)
end do
end if
end do
c
c The problem has been solved in the boundary, the vectors FI and DFI
c containing the values of the potential and its external derivatives
c along the boundary.
c
c Now it is possible to evaluate the function inside the domain. The
c integral eqn 2.5.10, Section 2.5, discretized in the form 3.5.20,
c Section 3.5.4, is applied to any point inside the domain previously
c specified.
c
do i = 1, npi
fint(i) = 0.0E+00
c
c Integrations along the boundary are performed from the point of
c interest. Numerical integration is now always used.
c
do j = 1, np
call hgnum ( xint(i),yint(i),x(j),y(j),x(j+1),y(j+1),
& a(1),a(2),b(1),b(2),igauss,xi,omeg)
if ( j .eq. np ) then
jj = 1
else
jj = j + 1
end if
fint(i) = fint(i) - a(1)*fi(j)-a(2)*fi(jj)+b(1)*dfi(j,2)+
& b(2)*dfi(jj,1)
end do
fint(i) = fint(i) / ( 2.0E+00 * 3.1415926E+00 )
end do
c
c The problem is finished and the results are going to be printed.
c
write(imp,1005) ienc
write(imp,1006) (i,fi(i),(dfi(i,j),j=1,2),i=1,np)
if ( 0 < npi ) then
write(imp,1007)
write(imp,1008) (xint(i),yint(i),fint(i),i=1,npi)
end if
c
c Terminate.
c
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BETIS:'
write ( *, '(a)' ) ' Normal end of execution.'
stop
c
c Wrong data
c
5 stop 5
11 stop 11
9988 write(*,*) 'Input File does not exist. '
stop 12
c
c Printing formats
c
1010 format(5x,a60,/,5x,30('**'))
1001 format(//,5x,'General constants :',/
@ /,5x,'Number of elements...............',i5,
@ /,5x,'Number of internal points........',i5,
@ /,5x,'Number of Gauss points...........',i5,
@ /,5x,'Node where the potential is null.',i5,//)
1002 format(5x,'Problem data:',
@ //,3x,'node',12x,'x coor',9x,'y coor',/)
1003 format(2x,i5,3x,2f15.4)
1102 format(/,3x,'node',5x,'code',6x,'potential',4x
@ ,'flux before',5x,'flux after',/)
1103 format(2x,i5,3x,i5,1x,3f15.5)
1005 format(///,5x,a60,/,5x,60('*'),//,5x,
@ 'SOLUTION ON THE BOUNDARY:',
@ //,5x,'boundary points'//6x,'node'
@ ,6x,'potential',4x,'flux before ',3x,'flux after',/)
1006 format(5x,i5,3f15.5)
1007 format(//,5X,'SOLUTION IN THE DOMAIN:',//,5x,
@ 'internal points',//,4x,'x coor',4x,'y coor',
@ 6x,'potential',/)
1008 format(2f10.3,f15.3)
end
subroutine hgnum ( xp, yp, x1, y1, x2, y2, a1, a2, b1, b2,
& igauss, xi, omeg )
c*********************************************************************72
c
cc HGNUM calculates integrals when the point is not in the element.
c
c Discussion:
c
c The routine calculates numerically the integrations along
c the elements of the boundary when the point where the integral
c equation is applied does not belong to the element, the general
c scheme represented by eqns 3.5.7, Section 3.5.2.1 being applied.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
c Parameters:
c
c Input, real XI(*), the natural coordinates of the Gauss
c points used in the numerical integration.
c
c Input, real OMEG(*), the weights associated with such points.
c
implicit none
real a1
real a2
real ax
real ay
real b1
real b2
real bx
real by
real dis
real g
real h
integer i
integer igauss
real omeg(50)
real r
real sig
real x1
real x2
real xc
real xi(50)
real xp
real y1
real y2
real yc
real yp
ax = ( x2 - x1 ) / 2.0E+00
ay = ( y2 - y1 ) / 2.0E+00
bx = ( x2 + x1 ) / 2.0E+00
by = ( y2 + y1 ) / 2.0E+00
c
c DIS represents the distance from the point where the integral
c equation is applied to the element along which the integration is
c being performed.
c
if ( ax .ne. 0.0E+00 ) then
dis = abs ( ( ay * xp / ax - yp + y1 - ay * x1 / ax )
& / sqrt ( ( ay / ax )**2 + 1.0E+00 ) )
else
dis = abs ( xp - x1 )
end if
sig = ( x1 - xp ) * ( y2 - yp ) - ( x2 - xp ) * ( y1 - yp )
if ( sig .lt. 0.0E+00 ) then
dis = - dis
end if
a1 = 0.0E+00
a2 = 0.0E+00
b1 = 0.0E+00
b2 = 0.0E+00
c
c The kernel of the integral is now evaluated at the Gauss points.
c
do i = 1, igauss
xc = ax * xi(i) + bx
yc = ay * xi(i) + by
c
c R represents the distance from the point where the integral
c equation is applied to the Gauss point under consideration.
c
r = sqrt ( ( xp - xc )**2 + ( yp - yc )**2 )
h = dis * omeg(i) * sqrt ( ax**2 + ay**2 ) / r**2
g = alog ( 1.0E+00 / r ) * omeg(i) * sqrt ( ax**2 + ay**2 )
c
c Equations 3.5.7.
c
a1 = a1 + ( xi(i) - 1.0E+00 ) * h / 2.0E+00
a2 = a2 - ( xi(i) + 1.0E+00 ) * h / 2.0E+00
b1 = b1 - ( xi(i) - 1.0E+00 ) * g / 2.0E+00
b2 = b2 + ( xi(i) + 1.0E+00 ) * g / 2.0E+00
end do
return
end
subroutine hgana ( x1, y1, x2, y2, a1, a2, b1, b2 )
c*********************************************************************72
c
cc HGANA calculates integrals along elements adjacent to the node.
c
c Discussion:
c
c Subroutine hgana is used to calculate analytically the integrations
c along the elements adjacent to the node where the integral equation
c is being applied.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
real a1
real a2
real b1
real b2
real dis
real x1
real x2
real y1
real y2
a1 = 0.0E+00
a2 = 0.0E+00
c
c The expressions 3.5.9 to 3.5.12, Section 3.5.2.2, are applied.
c
dis = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 )
b1 = dis * ( 1.5E+00 - alog ( dis ) ) / 2.0E+00
b2 = dis * ( 0.5E+00 - alog ( dis ) ) / 2.0E+00
return
end
subroutine pivo ( a, n, c, ipivo, npmax, imp, na )
c*********************************************************************72
c
cc PIVO applies Gauss elimination to solve the linear system.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
integer npmax
real a(npmax,npmax)
real aux
real c(npmax)
integer i
integer imp
integer ipivo(npmax)
integer j
integer k
integer l
real mx
integer n
integer na
integer pmx
do j = 1, n
mx = a(j,j)
pmx = j
do i = j + 1, n + na
if ( abs ( mx ) .lt. abs ( a(i,j) ) ) then
mx = a(i,j)
pmx = i
end if
end do
if ( abs ( mx ) .lt. 1.0E-06 ) then
write ( imp, '(a)' ) ' '
write ( imp, '(a)' ) 'PIVO - Fatal error!'
write ( imp, '(a)' ) ' The matrix is singular.'
stop 1111
else
if ( pmx .ne. j ) then
k = ipivo(pmx)
ipivo(pmx) = ipivo(j)
ipivo(j) = k
do l = j, n
aux = a(pmx,l)
a(pmx,l) = a(j,l)
a(j,l) = aux
end do
aux = c(pmx)
c(pmx) = c(j)
c(j) = aux
end if
do l = j + 1, n
a(j,l) = a(j,l) / a(j,j)
end do
c(j) = c(j) / a(j,j)
a(j,j) = 1.0E+00
do i = 1, n + na
if ( i .ne. j ) then
do l = j + 1, n
a(i,l) = a(i,l) - a(i,j) * a(j,l)
end do
c(i) = c(i) - a(i,j) * c(j)
a(i,j) = 0.0E+00
end if
end do
end if
end do
return
end
subroutine gauss_qn ( n, ep, om )
c*********************************************************************72
c
cc GAUSS_QN determines a Gauss quadrature rule.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
real*8 coef(50)
real ep(50)
real*8 epsilon(50)
integer expon
integer i
integer ii
integer ir0
integer ir1
integer k
integer n
real om(50)
real*8 omeg(50)
real*8 root(0:50,2)
real*8 tol
tol = 1.0E-14
if ( n .lt. 2 .or. 50 .lt. n ) then
stop 1111
end if
c
c Initialization.
c
ir1 = 1
ir0 = 2
root(0,1) = -1.0D+00
root(1,1) = 0.0D+00
root(2,1) = 1.0D+00
c
c Computation of polynomial coefficients and roots
c
do k = 2, n
do i = 0, n / 2
call coefic ( k, i, coef(i+1) )
end do
do i = 1, k
call roots ( k, coef, root(i-1,ir1), root(i,ir1),
& tol, root(i,ir0) )
end do
root(k+1,ir0) = 1.0D+00
root(0,ir0) = -1.0D+00
ii = ir0
ir0 = ir1
ir1 = ii
end do
c
c Computation of weights
c
do i = 1, n
epsilon(i) = root(i,ir1)
omeg(i) = 0.0D+00
do k = 0, n / 2
expon = n - 2 * k
if ( 1 .lt. expon ) then
omeg(i) = omeg(i) + coef(k+1)*expon*epsilon(i)**(expon-1)
else
omeg(i) = omeg(i) + coef(k+1) * expon
end if
end do
omeg(i) = omeg(i) / 2.0D+00**n
omeg(i) = 2.0D+00 / ( omeg(i)**2 * ( 1.0D+00 - epsilon(i)**2))
ep(i) = epsilon(i)
om(i) = omeg(i)
end do
return
end
subroutine coefic ( n, k, value )
c*********************************************************************72
c
cc COEFIC determines the coefficients of the polynomial defining the Gauss rule.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
integer i
integer inter
integer k
integer n
integer n1
integer n2
integer n3
real*8 sig
real*8 value
value = 1.0D+00
inter = k / 2
if ( inter * 2 - k .ne. 0 ) then
sig = -1.0D+00
else
sig = 1.0D+00
end if
n2 = n - k
n1 = 2 * n2
n3 = n - 2 * k
if ( n1 .eq. 0 ) then
value = 1.0D+00
else
do i = n2 + 1, n1
value = value * dble ( i )
end do
end if
if ( n3 .ne. 0 ) then
do i = 1, n3
value = value / dble ( i )
end do
end if
if ( k .ne. 0 ) then
do i = 1, k
value = value / dble ( i )
end do
end if
value = value * sig
return
end
subroutine evalua ( ncoef, coef, x, value )
c*********************************************************************72
c
cc EVALUA evaluates the polynomial defining the Gauss rule.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
integer ncoef
real*8 coef(ncoef)
integer expon
integer k
real*8 value
real*8 x
value = 0.0D+00
do k = 0, ncoef / 2
expon = ncoef - 2 * k
if ( expon .eq. 0 ) then
value = value + coef(k+1)
else
value = value + coef(k+1) * x**expon
end if
end do
return
end
subroutine roots ( ncoef, coef, x00, xff, error, xi )
c*********************************************************************72
c
cc ROOTS seeks roots of the polynomial defining the Gauss rule.
c
c Modified:
c
c 15 December 2007
c
c Author:
c
c Federico Paris, Jose Canas,
c
c Reference:
c
c Federico Paris, Jose Canas,
c Boundary Element Method: Fundamentals and Applications,
c Oxford, 1997,
c ISBN: 0-19-856543-7
c LC: TA347.B69.P34.
c
implicit none
integer ncoef
real*8 coef(ncoef)
real*8 errac
real*8 error
real*8 value0
real*8 valuef
real*8 valuei
real*8 x0
real*8 x00
real*8 xf
real*8 xff
real*8 xi
x0 = x00
xf = xff
call evalua ( ncoef, coef, x0, value0 )
call evalua ( ncoef, coef, xf, valuef )
errac = dabs ( xf - x0 )
do while ( error .lt. errac )
xi = ( xf + x0 ) / 2.0D+00
call evalua ( ncoef, coef, xi, valuei )
if ( dabs ( valuei ) .lt. error ) then
return
end if
if ( valuei * value0 .lt. 0.0D+00 ) then
xf = xi
valuef = valuei
else
x0 = xi
value0 = valuei
end if
errac = abs ( xf - x0 )
end do
return
end