subroutine DLELAPPER(N,N1,N2,tau,xm,ym,xb,yb,nx,ny,lg,
     & BCT,BCV,phi,dphi)

C     BVP governed by 2D Laplace's equation
C     Determining unknowns on the boundary by discontinuous linear elements
C     N: number of elements
C     tau*L is the distance of the two selected points from endpoints of each element
C     xm,ym: real arrays containing the two selected points on each of the elements
C     xb,yb: real arrays containing starting points of elements
C     nx,ny: real arrays containing normal vectors to elements
C     lg: real arrays containing lengths of elements
C     BCT: integer arrays containing type of boundary conditions
C     BCT(k)=0 if phi is specified on k-th element
C     BCV: real arrays containing values specified on the two points on each element
C     phi: real arrays containing values of phi on the two points on each element (returned)
C     dphi: real arrays containing values of normal derivative of phi on the two points of each element (returned)

      integer m,k,N,N1,N2,BCT(1000),k1

      double precision xm(1000),ym(1000),xb(1000),yb(1000),
     & nx(1000),ny(1000),lg(1000),BCV(1000),A(1000,1000),
     & B(1000),pi,PF1,PF2,delkm,delnkm,phi(1000),dphi(1000),F1,F2,
     & Z(1000),PF3,PF4,F3,F4,tau

      pi=4d0*datan(1d0)

      do 2 m=1,2*N+4*N1
	B(m)=0d0
	do 2 k=1,2*N+4*N1
	A(m,k)=0d0
2	continue                  !


      do 10 m=1,2*N
      k1=0
      B(m)=0d0
      do 5 k=1,N
      call DPF(xm(m),ym(m),xb(k),yb(k),nx(k),ny(k),lg(k),
     & PF1,PF2,PF3,PF4)
      F1=PF1/pi
      F2=PF2/pi
      F3=PF3/pi
      F4=PF4/pi
      if (k.eq.m) then
      delkm=1d0
      else
      delkm=0d0
      endif
      if (k.eq.(m-N)) then
      delnkm=1d0
      else
      delnkm=0d0
      endif
      if (BCT(k).eq.0) then
      A(m,k)=((1d0-tau)*F1-F3/lg(k))/(2d0*tau-1d0)
      A(m,N+k)=(-tau*F1+F3/lg(k))/(2d0*tau-1d0)
      B(m)=B(m)-((-(1d0-tau)*F2+F4/lg(k))
     & /(2d0*tau-1d0)-0.5d0*delkm)*BCV(k)
     & -((tau*F2-F4/lg(k))
     & /(2d0*tau-1d0)-0.5d0*delnkm)*BCV(N+k)
      else if (BCT(k).eq.1) then
      A(m,k)=(-(1d0-tau)*F2+F4/lg(k))/(2d0*tau-1d0)-0.5d0*delkm
      A(m,N+k)=(tau*F2-F4/lg(k))/(2d0*tau-1d0)-0.5d0*delnkm
      B(m)=B(m)-(((1d0-tau)*F1-F3/lg(k))*BCV(k)
     & +(-tau*F1+F3/lg(k))*BCV(N+k))/(2d0*tau-1d0)

      else
      k1=k1+1
      A(m,k)=(-(1d0-tau)*F2+F4/lg(k))/(2d0*tau-1d0)-0.5d0*delkm
      A(m,N+k)=(tau*F2-F4/lg(k))/(2d0*tau-1d0)-0.5d0*delnkm
      A(m,2*N+k1)=((1d0-tau)*F1-F3/lg(k))/(2d0*tau-1d0)
      A(m,2*N+2*N1+k1)=(-tau*F1+F3/lg(k))/(2d0*tau-1d0)

      endif
5     continue
10    continue

      do 12 k=1,N1
	A(2*N+k,k)=1d0
	A(2*N+N1+k,2*N1+N2+1-k)=-1d0
	A(2*N+N1+k,N+k)=1d0
	A(2*N+k,N+2*N1+N2+1-k)=-1d0

	A(2*N+2*N1+k,2*N+k)=1d0
	A(2*N+3*N1+k,2*N+2*N1+N2+1-k-N2)=1d0
	A(2*N+3*N1+k,2*N+2*N1+k)=1d0
	A(2*N+2*N1+k,2*N+2*N1+2*N1+N2+1-k-N2)=1d0

12    continue


      call solver(A,B,2*N+4*N1,1,Z)

      k1=0
      do 15 m=1,N
      if (BCT(m).eq.0) then
      phi(m)=BCV(m)
      phi(N+m)=BCV(N+m)
      dphi(m)=Z(m)
      dphi(N+m)=Z(N+m)
      else if (BCT(m).eq.1) then
      phi(m)=Z(m)
      phi(N+m)=Z(N+m)
      dphi(m)=BCV(m)
      dphi(N+m)=BCV(N+m)

      else
      k1=k1+1
      phi(m)=Z(m)
      phi(N+m)=Z(N+m)
      dphi(m)=Z(2*N+k1)
      dphi(N+m)=Z(2*N+2*N1+k1)

      endif
15    continue


      return
      end