c*******************************************************************
c     Purpose
c     =======
c     Implementation of the Natural Element Method (NEM) in 2D
c     (Isotropic, Linear Elasticity)
c
c     Author: N. Sukumar, TAM @ NU
c     Date: May 1997
c      
c*******************************************************************
c
      program nem2d
      implicit none
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/triangle.h"
#include "include/material.h"
#include "include/stiffness.h"
#include "include/fixity.h"
#include "include/traction.h"
#include "include/timing.h"
c
      real*8 constants(MAXCONST) ! material constants and problem type
c
      time_problem = etime(tarray)
      call input_data(constants,no_of_nodes,coord,fixity,uebc,
     &                node_boundary,no_of_triangles,tri_connect,
     &                tri_quad,tri_boundary,no_of_nbc_cells,
     &                node_nbc,tx,ty)
c
c     check if parameter settings and macro definitions are valid
c
      call check_parmac(no_of_nodes,no_of_triangles,no_of_nbc_cells)
c
c     set material constants and problem type (plane stress/strain)
c
      call material_constants(constants)
c
#ifdef PATCH_EBC
      call impose_ebc_patch()
#elif PLATE_HOLE
      call impose_ebc_hole()
#endif
c
c     compute circumcenter and square of the circumradius for all
c     the Delaunay triangles
c
      call compute_circum_centerrsq()
c
c     stiffness matrix assembly
c
c     initialize time for the assembly of the stiffness matrix
c
      time_stiff = dtime(tarray)
c
      call stiffness2d()
c
c     compute time for the assembly of the stiffness matrix
c
      time_stiff = dtime(tarray)
      write(*,*)'Time taken for K-matrix assembly:',time_stiff,' sec.'
c
c     external force assembly
c
c     initialize time for the assembly of the external force vector
c
      time_fext = dtime(tarray)
c
      call fext2d()
c
c     compute time for the assembly of the external force vector
c
      time_fext = dtime(tarray)
      write(*,*)'Time taken for ext. force assembly:',time_fext,' sec'
c
c     solution of the linear system: Au = b
c
c     initialize time for the solution of the linear system: Au = b
c
      time_solve = dtime(tarray)
c
      call solver2d()
c
c     compute time for the solution of the linear system: Au = b
c
      time_solve = dtime(tarray)
      write(*,*)'Time taken to solve the system Au=b:',time_solve,' sec'
c
c     output field variables at nodes and at Gauss points
c
      call output2d()
c
      time_problem = etime(tarray)
      write(*,*)'Time to run the problem: ',time_problem,' sec'
c
      stop
      end
c
c*******************************************************************
c     subroutine impose_ebc_patch()
c     Purpose
c     =======
c     Impose the essential boundary conditions for the displacment
c     patch test
c
c*******************************************************************
c
      subroutine impose_ebc_patch()
      implicit none
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/fixity.h"
c
      real*8 xI, yI, tol
      real*8 xmin, xmax, ymin, ymax
      data xmin /0.d0/
      data xmax /1.d0/
      data ymin /0.d0/
      data ymax /1.d0/
      integer i
c
      tol = 1.0d-6
      do 10 i = 1,no_of_nodes
         xI = coord(1,i) 
         yI = coord(2,i)
         if ((dabs(xI-xmin) .lt. tol) .or.
     &       (dabs(xI-xmax) .lt. tol) .or.
     &       (dabs(yI-ymin) .lt. tol) .or.
     &       (dabs(yI-ymax) .lt. tol)) then
            fixity(1,i) = 0
            fixity(2,i) = 0
            uebc(1,i) = xI
            uebc(2,i) = yI
         endif
10    continue
c
      return
      end
c
c*******************************************************************
c     subroutine impose_ebc_hole()
c     Purpose
c     =======
c     Impose the symmetry EBCs for a plate with a hole
c
c*******************************************************************
c
      subroutine impose_ebc_hole()
      implicit none
c
c     include files
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/fixity.h"
c
      real*8 xI, yI, tol
      integer i
c
      tol = 1.0d-6
      do 10 i = 1,no_of_nodes
         xI = coord(1,i) 
         yI = coord(2,i)
         if (dabs(xI) .lt. tol) then
            fixity(1,i) = 0
            uebc(1,i) = 0.d0
         endif
         if (dabs(yI) .lt. tol) then
            fixity(2,i) = 0
            uebc(2,i) = 0.d0
         endif
10    continue        
c
      return
      end
