c**************************************************************
c     subroutine fext2d()
c     Purpose
c     =======
c     Subroutine to assemble the external force vector for the
c     Natural Element Method (NEM)
c
c**************************************************************
c
      subroutine fext2d()
      implicit none
c
c     include files for the external force assembly (NEM program)
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/gausspt_x.h"
#include "include/fixity.h"
#include "include/traction.h"
#include "include/penalty.h"
#include "include/fext.h"
c
c     declarations for the external force assembly
c
      real*8 xynbc(NSD,NODES_PER_NBCCELL)
      real*8 N_mat(NSD,NODES_PER_NBCCELL)
      real*8 sj, da
      real*8 psi(15), weight(15)
      integer nintnbc, kx
      integer nbc_cell_no, index
      integer nodes(NODES_PER_NBCCELL)
      real*8 traction(NSD,NODES_PER_NBCCELL), txsum, tysum
      real*8 row, rows(NSD), fext(NSD)
c
      integer i, j
c
c     initialize external force vector to zero
c
      do 05 i = 1,NSD*no_of_nodes
	 b(i) = 0.d0
05    continue
c
c     compute the external force vector b [b is initialized to
c     zero in nem.F]
c
      do 10 nbc_cell_no = 1,no_of_nbc_cells
	 index = (NODES_PER_NBCCELL)*nbc_cell_no - NODES_PER_NBCCELL
	 do 20 i = 1,NODES_PER_NBCCELL
            nodes(i) = node_nbc(index+i)
	    traction(1,i) = tx(index+i)
	    traction(2,i) = ty(index+i)
20       continue
c
	 nintnbc = 3   ! fixing the quadrature for the integration
c                      ! along the natural boundary \Gamma_t
c
c     get psi(nintnbc) and weight(nintnbc) -- natural coordinates
c     (xi,eta) and associated weights
c
         call gqb(psi,weight,nintnbc)
c
	 do 30 i = 1,NODES_PER_NBCCELL
	    xynbc(1,i) = coord(1,nodes(i))
	    xynbc(2,i) = coord(2,nodes(i))
30       continue
c
         n = 2
         ngbrs_list(1) = nodes(1)
         ngbrs_list(2) = nodes(2)
c
         do 40 kx = 1,nintnbc
c
c     compute the Jacobian and linear (FE) shape functions
c     N_mat: N_mat(i,j) [i=1,2; j=1,2]; sj: Jacobian
c
            N_mat(1,1) = (1.0 - psi(kx))/2.0
            N_mat(1,2) = (1.0 + psi(kx))/2.0
            xbar = N_mat(1,1)*xynbc(1,1) + N_mat(1,2)*xynbc(1,2)
            ybar = N_mat(1,1)*xynbc(2,1) + N_mat(1,2)*xynbc(2,2)
	    sj = 0.5d0*(dsqrt((xynbc(1,2) - xynbc(1,1))**2 +
     &      	            (xynbc(2,2) - xynbc(2,1))**2))
            da = sj*weight(kx)
c
c     evaluate the components of the traction vector at the Gauss 
c     point: t = \sum_{i=1}^{NODES_PER_NBCCELL} N_i t_i 
c
c     tractions [natural BCs] for plate with a hole problem
c
#ifdef PLATE_HOLE
            call compute_nbc_hole(xbar,ybar,txsum,tysum)
#else
c
            txsum = 0.0d0
            tysum = 0.0d0
            do 50 j = 1,NODES_PER_NBCCELL
               txsum = txsum + N_mat(1,j)*traction(1,j)
               tysum = tysum + N_mat(1,j)*traction(2,j)
50          continue
#endif
c
c     natural boundary conditions: FE shape functions
c
            phi(1,1) = N_mat(1,1)
            phi(1,2) = N_mat(1,2)
c
            do 60 i = 1,n
               rows(1) = NSD*ngbrs_list(i) - 1
 	       rows(2) = NSD*ngbrs_list(i)
               fext(1) = phi(1,i)*txsum*da
               fext(2) = phi(1,i)*tysum*da
c
c     add the nodal external force vector to the RHS vector b
c
               b(rows(1)) = b(rows(1)) + fext(1)
               b(rows(2)) = b(rows(2)) + fext(2)
60          continue          
40       continue
10    continue
c
c     impose the essential BCs in the RHS vector: if u_i is known,
c     then set the ith entry of the external force vector to a
c     a very large +ve real number (penalty parameter) multiplied
c     by the imposed essential BC
c
      do 70 i = 1,no_of_nodes 
         if (fixity(1,i) .eq. 0) then
            row = NSD*i - 1
            b(row) = uebc(1,i)*PENALTY_PARAMETER
         endif
         if (fixity(2,i) .eq. 0) then
            row = NSD*i
            b(row) = uebc(2,i)*PENALTY_PARAMETER
         endif
70    continue
c
      return
      end
c
c*******************************************************************
c     subroutine compute_nbc_hole(x,y,tx,ty)
c     Purpose
c     =======
c     Compute the exact tractions (plate with a hole - plate has
c     infinite extent in x- and y-directions) on the boundary. Unit
c     traction (at infinity) in the x-direction is imposed
c
c*******************************************************************
c
      subroutine compute_nbc_hole(x,y,tx,ty)
      implicit none
#include "include/plate_hole.h"
c
      real*8 x, y, tx, ty
      real*8 tol, r, theta, aovr, sigx, sigy, sigxy
c
      tol = 1.0d-3
      r = dsqrt(x**2 + y**2)
      theta = datan2(y,x)
      aovr = a/r
      sigx = 1. - aovr**2 * (1.5 * dcos(2.*theta)
     &       + dcos(4.*theta)) + 1.5*aovr**4 * dcos(4.*theta)
      sigy = -aovr**2 * (0.5 * dcos(2.*theta) - dcos(4.*theta))
     &       - 1.5*aovr**4 * dcos(4.*theta)
      sigxy= -aovr**2 * (0.5 * dsin(2.*theta) + dsin(4.*theta))
     &       + 1.5*aovr**4 * dsin(4.*theta)
c
      if (dabs(x-5.) .lt. tol) then
         tx = sigx
         ty = sigxy
      endif
      if (dabs(y-5.) .lt. tol) then
         tx = sigxy
         ty = sigy
      endif
c
      return
      end
