PROGRAM num_int
!-------------------------------------------------------------------------
! Numerical integration of finite elements
!-------------------------------------------------------------------------
 USE main; IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::i,j,k,ndim=3,ndof,nip,nod,nodof=3,nprops=2,nst=6
 REAL(iwp)::det,zero=0.0_iwp,start_time,end_time
 CHARACTER(len=15)::element
!-----------------------dynamic arrays------------------------------------
 REAL(iwp),ALLOCATABLE::bee(:,:),coord(:,:),dee(:,:),der(:,:),deriv(:,:), &
   jac(:,:),km(:,:),points(:,:),prop(:,:),weights(:)  
!-----------------------input and initialisation--------------------------
 OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
 READ(10,*)element,nod,nip
 IF(element=='triangle'.OR.element=='quadrilateral')THEN
   ndim=2; nodof=2; nst=3
 END IF
 ndof=nod*nodof
 ALLOCATE(points(nip,ndim),dee(nst,nst),coord(nod,ndim),jac(ndim,ndim),   &
   weights(nip),der(ndim,nod),deriv(ndim,nod),bee(nst,ndof),km(ndof,ndof),&
   prop(nprops,1))
 READ(10,*)prop(:,1)
 WRITE(11,'("   NUMERICAL INTEGRATION")')
 WRITE(11,'("            OF")')
 WRITE(11,'(" ELEMENT STIFFNESS MATRICES",/)')
 WRITE(11,'(i2,"-node",2x,a)')nod,element
 WRITE(11,'(i2,"-integrating point(s)")')nip
 WRITE(11,'(/," Nodal (x,y[,z]) coordinates")')
 DO i=1,nod
   READ(10,*)coord(i,:)
   WRITE(11,'(3f12.4)')coord(i,:)
 END DO
 WRITE(11,'(/," Youngs modulus and Poissons ratio")')
 WRITE(11,'(2f12.4/)')prop(1:2,1)
 WRITE(11,'(/," Stiffness matrix")')
!-----------------------element stiffness integration and assembly--------
 CALL sample(element,points,weights)
 CALL deemat(dee,prop(1,1),prop(2,1))
 Call cpu_time(start_time)
! DO j=1,1000000
 km=zero
 int_pts_1: DO i=1,nip
   CALL shape_der(der,points,i)
   jac=MATMUL(der,coord); det=determinant(jac); CALL invert(jac)
   deriv=MATMUL(jac,der); CALL beemat(bee,deriv)
   km=km+MATMUL(MATMUL(transpose(bee),dee),bee)*det*weights(i)
 END DO int_pts_1
! END DO
 DO i=1,nod*ndim
   DO k=1,i
     WRITE(11,'(2i3,f20.4)')i,k,km(i,k)
   END DO
 END DO 
 Call cpu_time(end_time)
 WRITE(11,'(/,A,f12.4)')"time taken=",end_time-start_time
STOP
END PROGRAM num_int
