c**************************************************************
c     subroutine stiffness2d()
c     Purpose
c     =======
c     Subroutine to assemble the stiffness matrix for the
c     Natural Element Method (NEM)
c
c**************************************************************
c
      subroutine stiffness2d()
      implicit none 
c
c     include files for the stiffness assembly (NEM program)
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"  ! declaration for n-n arrays
#include "include/triangle.h"
#include "include/gausspt_x.h"
#include "include/gausspt_ref.h"
#include "include/quadrature.h"
#include "include/material.h"
#include "include/stiffness.h"
#include "include/fext.h"
#include "include/fixity.h"
#include "include/penalty.h"
c
c     determinant of the Jacobian and the area
c
      real*8 detj, dv
c
c     nodes and nodal coordinates of the vertices of a triangle
c
      integer nd1, nd2, nd3
      real*8 vertices(NSD,NODES_PER_TRIANGLE)
c
c     declarations for the nodal-stiffness matrix
c
      integer row, col
      real*8 val1, val2, val3, val4
c
c     local integer counters
c
      integer i, j, ngqp, tri_no
c
      integer count
      count = 0
c
c     initialize stiffness matrix to zero
c
      do 10 i = 1,NSD*no_of_nodes
         do 20 j = 1,NSD*no_of_nodes
            A(j,i) = 0.d0
20       continue
10    continue
c
c     assembly of the stiffness matrix
c
c     loop over all the triangles
c
      do 30 tri_no = 1,no_of_triangles
c
c     get pgauss(3,nint) and wgauss(nint) - natural coordinates and
c     associated weights
c
         nint = tri_quad(tri_no)
         call gaussq_triangle(nint,pgauss,wgauss)
c
c     get nodal coordinates of the nodes at the vertex of triangle
c     "tri_no"
c
         nd1 = tri_connect(1,tri_no)
         nd2 = tri_connect(2,tri_no)
         nd3 = tri_connect(3,tri_no)
         vertices(1,1) = coord(1,nd1)
         vertices(2,1) = coord(2,nd1)
         vertices(1,2) = coord(1,nd2)
         vertices(2,2) = coord(2,nd2)
         vertices(1,3) = coord(1,nd3)
         vertices(2,3) = coord(2,nd3)
         open(11,file='tricount.dat')
         write(11,*)'Processing triangle no',tri_no
         close(11)
c
         do 40 ngqp = 1,nint
            xi1 = pgauss(1,ngqp)
            xi2 = pgauss(2,ngqp)
            xi3 = pgauss(3,ngqp)
c
c     get Jacobian detj and the global coordinates (xbar,ybar) 
c     associated with the quadrature point xi1,xi2,xi3
c
            call compute_xgJ_triangle(vertices,detj)
            dv = 0.5*detj*wgauss(ngqp)
c
#ifdef FE
c
c     compute FE shape functions and derivatives
c
            call feshape2d_tri(tri_no,n,ngbrs_list,phi)
#else
c
c     compute number of natural neighbors for the Gauss point,
c     the neighbor list in ngbrs_list(), the NEM shape functions
c     (n-n coordinates) phi_I and its derivatives phi_I,j
c
            call nemshape2d(n,ngbrs_list,phi)
#endif
c
            do 50 I = 1,n
               do 60 J = 1,n
                  if (ngbrs_list(I) .le. ngbrs_list(J)) then
c
c     insert nodal stiffness matrix into the global stiffness 
c     matrix A(.,.)
c
                     val1 = phi(2,I)*phi(2,J)
                     val2 = phi(3,I)*phi(3,J)
                     val3 = phi(2,I)*phi(3,J)
                     val4 = phi(3,I)*phi(2,J)
                     row = NSD*ngbrs_list(I)
                     col = NSD*ngbrs_list(J)
                     A(row-1,col-1) = A(row-1,col-1) +
     &                                (val1*d1 + val2*d3)*dv
                     A(row,col) = A(row,col) + 
     &                            (val1*d3 + val2*d1)*dv
                     A(row-1,col) = A(row-1,col) +
     &                              (val3*d2 + val4*d3)*dv
                     if (ngbrs_list(I) .ne. ngbrs_list(J)) then
                        A(row,col-1) = A(row,col-1) +
     &                                 (val3*d3 + val4*d2)*dv
                     endif
                  else
c
c     skip the lower triangle of the stiffness matrix
c
                     continue
                  endif
60             continue
50          continue
40       continue
30    continue
c
c     impose the essential BCs by inserting a penalty parameter on
c     the diagonal: if u_i is known, then a very large positive 
c     real number (penalty parameter) is inserted on the diagonal
c
      open(60,file='EBCs.dat')
      write(60,*)'INFO ON EBCs'
      write(60,*)'============'
      write(60,*)
      write(60,*)'Row #   Node #    Coord_x   Coord_y   Displacement' 
      write(60,*)
      do 70 i = 1,no_of_nodes
         if (fixity(1,i) .eq. 0) then
            row = NSD*i - 1
            A(row,row) = PENALTY_PARAMETER ! row = col for the diag.
            write(60,199)row,i,coord(1,i),coord(2,i),uebc(1,i)
         endif
         if (fixity(2,i) .eq. 0) then
            row = NSD*i
            A(row,row) = PENALTY_PARAMETER ! row = col for the diag.
            write(60,199)row,i,coord(1,i),coord(2,i),uebc(2,i)
         endif
70    continue
199   format(1x,I5,5x,I5,5x,F15.10,5x,F15.10,5x,F15.10)
      close(60)
c
      return
      end
