c**************************************************************
c     subroutine nemshape2d(n,ngbrs_list,phi)
c     Purpose
c     =======
c     Subroutine to compute the number of natural neighbors for
c     a Gauss/sampling point, the nodal neighbor list, the
c     shape functions phi_I [I=1,n] (natural neighbor coords.),
c     and the derivatives of the shape functions phi_I,j [I=1,n
c     and j=1,2] (derivatives of the natural neighbor coords.)
c     for the Natural Element Method (NEM)
c
c     Author: N. Sukumar, TAM @NU
c     Date:   May 1997
c**************************************************************
c
      subroutine nemshape2d(n,ngbrs_list,phi)
      implicit none
c
c     include files for the n-n coordinates (NEM shape functions)
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/gausspt_x.h"
c
c     pointer array for the position of the neighbors in array
c     ngbrs_list()
c
      integer ngbrptr_list(no_of_nodes)
c
c     vertices (nodes) of a triangle
c
      integer nd(NODES_PER_TRIANGLE)
c
      real*8 a1, a2, b1, b2, p1, p1x, p1y, p2, p2x, p2y
c
c     circumcenter of a triangle
c      
      real*8 v1, v2
c
c     circumcenter coordinates c and their derivatives c,i for the
c     circum-triangles
c
      real*8 c1(NSD+1,NODES_PER_TRIANGLE)
      real*8 c2(NSD+1,NODES_PER_TRIANGLE)
c
      real*8 dar2(NSD+1), totarea2(NSD+1)
c
c     integer counters
c
      integer tri_no, tri_global_no, i, j, k
c
      call find_neighbors(n,ngbrs_list,ngbrptr_list,ntri,
     &                    ngbrs_tri)
c
      do 05 i = 1,n
         phi(1,i) = 0.d0
         phi(2,i) = 0.d0
         phi(3,i) = 0.d0
05    continue  
c
      totarea2(1) = 0.d0
      totarea2(2) = 0.d0
      totarea2(3) = 0.d0
c
      do 10 tri_no = 1,ntri
         tri_global_no = ngbrs_tri(tri_no)
         v1 = tri_ccenter(1,tri_global_no)
         v2 = tri_ccenter(2,tri_global_no)
         do 20 i = 1,NODES_PER_TRIANGLE
            if (i .eq. 1) then
               a1 = coord(1,tri_connect(2,tri_global_no))
               a2 = coord(2,tri_connect(2,tri_global_no))
               b1 = coord(1,tri_connect(3,tri_global_no))
               b2 = coord(2,tri_connect(3,tri_global_no))
            elseif (i .eq. 2) then
               a1 = coord(1,tri_connect(3,tri_global_no))
               a2 = coord(2,tri_connect(3,tri_global_no))
               b1 = coord(1,tri_connect(1,tri_global_no))
               b2 = coord(2,tri_connect(1,tri_global_no)) 
            else
               a1 = coord(1,tri_connect(1,tri_global_no))
               a2 = coord(2,tri_connect(1,tri_global_no))
               b1 = coord(1,tri_connect(2,tri_global_no))
               b2 = coord(2,tri_connect(2,tri_global_no))
            endif
            call circumcenter(a1,a2,b1,b2,xbar,ybar,p1,p2,p1x,p1y,
     &                        p2x,p2y)
            c1(1,i) = p1
            c1(2,i) = p1x
            c1(3,i) = p1y
            c2(1,i) = p2
            c2(2,i) = p2x
            c2(3,i) = p2y
20       continue 
         do 30 j = 1,NODES_PER_TRIANGLE
            nd(j) = tri_connect(j,tri_global_no)
            if (j .eq. 1) then
               call darea2(c1(1,2),c2(1,2),c1(1,3),c2(1,3),v1,v2,dar2)
            elseif (j .eq. 2) then
               call darea2(c1(1,3),c2(1,3),c1(1,1),c2(1,1),v1,v2,dar2)
            else
               call darea2(c1(1,1),c2(1,1),c1(1,2),c2(1,2),v1,v2,dar2)
            endif
c
            i = ngbrptr_list(nd(j))
            do 40 k = 1,NSD+1
               phi(k,i) = phi(k,i) + dar2(k)
               totarea2(k) = totarea2(k) + dar2(k)
40          continue
30       continue
10    continue
c
      if (totarea2(1) .gt. 0.d0) then
         do 50 i = 1,n
            phi(1,i) = phi(1,i)/totarea2(1)
            phi(2,i) = (phi(2,i) - phi(1,i)*totarea2(2))/totarea2(1)
            phi(3,i) = (phi(3,i) - phi(1,i)*totarea2(3))/totarea2(1)
50       continue
      else
         write(6,*)'Error: Shape functions are -ve/zero:',xbar,ybar
         write(6,*)'Neighbors are',(ngbrs_list(i),i=1,n)
         write(6,*)'Circum-triangles are',(ngbrs_tri(i),i=1,ntri)
         do i = 1,ntri
            write(6,*)'Connectivity for tri',i,' :',
     &                (tri_connect(j,ngbrs_tri(i)),j=1,3)
         enddo
         write(6,*)'xbar, ybar, total area:',xbar,ybar,totarea2(1)
         stop
      endif
c
c     check constant and linear consistency
c
      call consistency(n,ngbrs_list,phi)
c
      return
      end
c
c**************************************************************
c     subroutine find_neighbors(n,ngbrs_list,ngbrptr_list,ntri,
c    &                          ngbrs_tri)
c     Purpose
c     =======
c     Subroutine to determine the natural neighbors and all the
c     associated triangles which have vertices that are natural
c     neighbors of the point X(xbar,ybar)
c
c**************************************************************
c
      subroutine find_neighbors(n,ngbrs_list,ngbrptr_list,ntri,
     &                          ngbrs_tri)
      implicit none  
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/gausspt_x.h"  ! Gauss point (xbar,ybar) 
c
      integer ngbrptr_list(*)
c
      real*8 rsq_xg
c
c     vertices (nodes) of a triangle
c
      integer nd1, nd2, nd3
c
c     integer counters
c
      integer i
c
      do 10 i = 1,no_of_nodes
         ngbrptr_list(i) = 0
10    continue
c  
      ntri = 0
      do 20 i = 1,no_of_triangles
c
         rsq_xg = (tri_ccenter(1,i) - xbar)*(tri_ccenter(1,i) - xbar)
     &            + (tri_ccenter(2,i) - ybar)*(tri_ccenter(2,i) - ybar)
         if (rsq_xg .le. tri_sqradius(i)) then
c
c     (xbar,ybar) is inside the circumcircle about triangle "i"
c
            ntri = ntri + 1
            ngbrs_tri(ntri) = i
c
            nd1 = tri_connect(1,i)
            nd2 = tri_connect(2,i)
            nd3 = tri_connect(3,i)
            ngbrptr_list(nd1) = 1
            ngbrptr_list(nd2) = 1
            ngbrptr_list(nd3) = 1
         endif
20    continue
c
      n = 0
      do 30 i = 1,no_of_nodes
         if (ngbrptr_list(i) .eq. 1) then
            n = n + 1
            ngbrs_list(n) = i
            ngbrptr_list(i) = n
         endif
30    continue
c
      if (n .gt. MAXNGBRS) then
         write(6,*)'Error: Increase MAXNGBRS >=',n
         stop
      endif 
      if (ntri .gt. MAXTRINGBRS) then
         write(6,*)'Error: Increase MAXTRINGBRS >=',ntri
         stop
      endif
c
      return
      end
c
c**************************************************************
c     subroutine circumcenter(a1,a2,b1,b2,x,y,v1,v2,v1x,v1y,
c    &                        v2x,v2y)
c     Purpose
c     =======
c     Subroutine to determine the circumcenter (v1,v2) and the
c     its derivatives v1,i and v2,i (i = 1,2) for the triangle
c     (p_1,p_2,X), with vertices p_1(a1,a2), p_2(b1,b2) and 
c     X(x,y)
c
c**************************************************************
c
      subroutine circumcenter(a1,a2,b1,b2,x,y,v1,v2,v1x,v1y,
     &                        v2x,v2y)
      implicit none
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/triangle.h"
c
      real*8 TOL
      parameter (TOL = 1.d-18)
      real*8 a1, a2, b1, b2, x, y  ! coordinates of the vertices
      real*8 v1, v2, v1x, v1y, v2x, v2y 
      real*8 term1, term2, term3
      real*8 D, Dx, Dy ! D is the determinant (Cramer's rule)
c
      D = (a1 - x)*(b2 - y) - (b1 - x)*(a2 - y)
      if (dabs(D) .lt. TOL) then
         write(6,*)'Error:Three points of the triangle are collinear'
         write(6,*)a1,a2,b1,b2,x,y
         stop
      endif
      term1 = 0.5*(((a1 - x)*(a1 + x) + (a2 - y)*(a2 + y))*(b2 - y)
     &        - ((b1 - x)*(b1 + x) + (b2 - y)*(b2 + y))*(a2 - y))
      v1 = term1/D
      term2 = 0.5*(((b1 - x)*(b1 + x) + (b2 - y)*(b2 + y))*(a1 - x)
     &        - ((a1 - x)*(a1 + x) + (a2 - y)*(a2 + y))*(b1 - x))
      v2 = term2/D
c
      Dx = a2 - b2
      Dy = b1 - a1
      term3 = 0.5d0*((b1 + a1)*(b1 - a1) + (b2 + a2)*(b2 - a2))
      v1x = (x - v1)*Dx/D
      v1y = (term3 + y*Dx - v1*Dy)/D
      v2x = (-term3 + x*Dy - v2*Dx)/D
      v2y = (y - v2)*Dy/D
c
      return
      end
c
c**************************************************************
c     subroutine compute_circum_centerrsq()
c     Purpose
c     =======
c     Subroutine to compute the circumcenter and square of the
c     circum-radius for all the Delaunay triangles
c
c**************************************************************
c
      subroutine compute_circum_centerrsq()
      implicit none
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/triangle.h"
c
      real*8 a1, a2, b1, b2, c1, c2  ! coordinates of the vertices
      real*8 D  ! D is the determinant (Cramer's rule)
c
c     vertices (nodes) of a triangle
c
      integer nd1, nd2, nd3
c
c     integer counter
c
      integer tri_no
c
c 
      do 10 tri_no = 1,no_of_triangles
         nd1 = tri_connect(1,tri_no)
         nd2 = tri_connect(2,tri_no)
         nd3 = tri_connect(3,tri_no)
         a1 = coord(1,nd1)
         a2 = coord(2,nd1)
         b1 = coord(1,nd2)
         b2 = coord(2,nd2)
         c1 = coord(1,nd3)
         c2 = coord(2,nd3)
c
         D = (a1 - c1)*(b2 - c2) - (b1 - c1)*(a2 - c2)
c
         tri_ccenter(1,tri_no) = 0.5*(((a1 - c1)*(a1 + c1) + 
     &                                (a2 - c2)*(a2 + c2))*(b2 - c2)
     &                                - ((b1 - c1)*(b1 + c1) +
     &                                (b2 - c2)*(b2 + c2))*(a2 - c2))/D
         tri_ccenter(2,tri_no) = 0.5*(((b1 - c1)*(b1 + c1) + 
     &                                (b2 - c2)*(b2 + c2))*(a1 - c1)
     &                                - ((a1 - c1)*(a1 + c1) + 
     &                                (a2 - c2)*(a2 + c2))*(b1 - c1))/D
         tri_sqradius(tri_no) = (c1 - tri_ccenter(1,tri_no))
     &                          *(c1 - tri_ccenter(1,tri_no)) +
     &                          (c2 - tri_ccenter(2,tri_no))
     &                          *(c2 - tri_ccenter(2,tri_no))
10    continue  
c
      return
      end
c
c**************************************************************
c     subroutine darea2(a1,a2,b1,b2,c1,c2,dar2)
c     Purpose
c     =======
c     Compute twice the (signed) area and its derivatives for
c     the triangle with vertices p_1(a1,a2), p_2(b1,b2) and 
c     p_3(c1,c2)
c
c**************************************************************
c
      subroutine darea2(a1,a2,b1,b2,c1,c2,dar2)
      implicit none
c
      real*8 a1(*), a2(*), b1(*), b2(*), c1, c2
      real*8 dar2(*)
c
c     a1(), a2(), b1(), and b2() are functions of x
c
      dar2(1) = (a1(1) - c1)*(b2(1) - c2) 
     &          - (b1(1) - c1)*(a2(1) - c2)
      dar2(2) = (a1(1) - c1)*b2(2) + (b2(1) - c2)*a1(2) 
     &          - (b1(1) - c1)*a2(2) - (a2(1) - c2)*b1(2) 
      dar2(3) = (a1(1) - c1)*b2(3) + (b2(1) - c2)*a1(3) 
     &          - (b1(1) - c1)*a2(3) - (a2(1) - c2)*b1(3)
c
      return
      end
c
c**************************************************************
c     subroutine boundary_triangles()
c     Purpose
c     =======
c     Subroutine to determine all the triangles t_i that have
c     least one node on the boundary of the convex hull
c
c**************************************************************
c
      subroutine boundary_triangles()
      implicit none  
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/triangle.h"
c
      integer nd(NODES_PER_TRIANGLE)
c
      logical boundary
c
c     integer counters
c
      integer i, j
c
      do 10 i = 1,no_of_triangles
         j = 1
         boundary = .false.
         do 20 while ((j .le. 3) .and. (.not. boundary))
            nd(j) = tri_connect(j,i)
            if (node_boundary(nd(j)) .eq. 1) then
               boundary = .true.
               tri_boundary(i) = 1
            endif
            j = j + 1
20       enddo        
10    continue  
c
      return
      end
c
c**************************************************************
c     subroutine consistency(n,ngbrs_list,phi)
c     Purpose
c     =======
c     Check constant and linear consistency
c
c**************************************************************
c
      subroutine consistency(n,ngbrs_list,phi)
      implicit none
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/gausspt_x.h"
c
      integer i
      real*8 sum, sumx, sumy, tol, sumphix, sumphiy, sumderx, sumdery
      real*8 sumder_x, sumder_y
      integer consistency_flag
      common /c_flag/ consistency_flag
c
      consistency_flag = 0
      tol = 1.0e-10
      sum = 0.d0
      sumx = 0.d0
      sumy = 0.d0
      sumphix = 0.d0
      sumphiy = 0.d0
      sumderx = 0.d0
      sumdery = 0.d0
      sumder_x = 0.d0
      sumder_y = 0.d0
      do 10 i = 1,n
         sum = sum + phi(1,i)
         sumx = sumx + phi(1,i)*coord(1,ngbrs_list(i))
         sumy = sumy + phi(1,i)*coord(2,ngbrs_list(i))
         sumphix = sumphix + phi(2,i)
         sumphiy = sumphiy + phi(3,i)
         sumderx = sumderx + phi(2,i)*coord(1,ngbrs_list(i))
         sumdery = sumdery + phi(3,i)*coord(2,ngbrs_list(i))
         sumder_x = sumder_x + phi(3,i)*coord(1,ngbrs_list(i))
         sumder_y = sumder_y + phi(2,i)*coord(2,ngbrs_list(i))
10    continue
c
      if ((dabs(1.d0-sum) .gt. tol) .or. (dabs(xbar-sumx) .gt. tol) 
     &    .or. (dabs(ybar-sumy) .gt. tol)
     &    .or. (dabs(sumphix) .gt. tol)
     &    .or. (dabs(sumphiy) .gt. tol)
     &    .or. (dabs(sumdery-1.d0) .gt. tol)
     &    .or. (dabs(sumdery-1.d0) .gt. tol)
     &    .or. (dabs(sumder_x) .gt. tol)
     &    .or. (dabs(sumder_y) .gt. tol)) then
      write(6,*)'Consistency error at',xbar,ybar,dsqrt(xbar**2+ybar**2)
         write(6,*)sum,sumx-xbar,sumy-ybar
         write(6,*)sumphix,sumphiy
         write(6,*)sumderx-1.d0,sumdery-1.d0,sumder_x,sumder_y
         write(6,*)
         consistency_flag = 1
      endif
c
      return
      end
