!     Program to test integration of derivatives of Tetra1
!     ----------------------------------------------------

!     From the values of some functions and their first derivatives (v) at the 
!     nodes of a random tetrahedron (coor), we compute the integrals by using 
!     the shape functions (phi) defined at a set of integration points (point) 
!     and the related weights (weight).

      IMPLICIT INTEGER*4 (i-n), REAL*8 (a-h,o-z)

      DIMENSION coor(3,4),point(3,31),phi(4,16,31),weight(31),v(4)
      REAL*4 xx

      OPEN (UNIT=10,FILE='resu2.txt',STATUS='UNKNOWN')

!     random tetra

      iran = 1.
      DO i = 1,3
         DO j = 1,4
            xx = RAND (iran);  iran = 0
            coor(i,j) = xx*10.d0
         ENDDO
      ENDDO
      x1 = coor(1,1);  y1 = coor(2,1);  z1 = coor(3,1)
      x2 = coor(1,2);  y2 = coor(2,2);  z2 = coor(3,2)
      x3 = coor(1,3);  y3 = coor(2,3);  z3 = coor(3,3)
      x4 = coor(1,4);  y4 = coor(2,4);  z4 = coor(3,4)

      WRITE (10,4)
    4 FORMAT (' Hermite (Délèze) tetrahedron interpolation with 16 parameters'/            &
              ' Integration of cubic terms on a random tetrahedron.'/                      &
              ' Check integration of the derivatives.'/)

      WRITE (10,'('' x1,y1,z1 = '',1P3D15.6)') x1,y1,z1
      WRITE (10,'('' x2,y2,z2 = '',1P3D15.6)') x2,y2,z2
      WRITE (10,'('' x3,y3,z3 = '',1P3D15.6)') x3,y3,z3
      WRITE (10,'('' x4,y4,z4 = '',1P3D15.6)') x4,y4,z4

!     shape functions

      npi = 5
      ivers = 1
   10 WRITE (10,12) npi
   12 FORMAT (/' Nb. of integration points = ',i3)
      CALL FElib_Tetra_Shape (coor,npi,point,ivers,phi,weight)

!     check integration of u(x,y,z) = x

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(1,j);  v(2) = 1.d0;  v(3) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ &
                       phi(1,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      WRITE (10,'('' Int (x)               = '',1PD15.8)') sum

!     check integration of u(x,y,z) = 0.5d(x**2)/dx

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(1,j)**2;  v(2) = 2*coor(1,j);  v(3) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(2,j1,i)*v(1)+phi(2,j1+1,i)*v(2)+phi(2,j1+2,i)*v(3)+ &
                       phi(2,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      sum = sum/2.d0
      WRITE (10,'('' Int (0.5d(x**2)/dx)   = '',1PD15.8)') sum

!     check integration of u(x,y,z) = x**2

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(1,j)**2;  v(2) = 2*coor(1,j);  v(3) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ &
                       phi(1,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      WRITE (10,'('' Int (x**2))           = '',1PD15.8)') sum

!     check integration of u(x,y,z) = 0.33*d(x**3)/dx

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(1,j)**3;  v(2) = 3*coor(1,j)**2;  v(3) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(2,j1,i)*v(1)+phi(2,j1+1,i)*v(2)+phi(2,j1+2,i)*v(3)+ &
                       phi(2,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      sum = sum/3.d0
      WRITE (10,'('' Int (0.33*d(x**3)/dx) = '',1PD15.8)') sum

!     check integration of u(x,y,z) = y**2

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(2,j)**2;  v(3) = 2*coor(2,j);  v(2) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ &
                       phi(1,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      WRITE (10,'('' Int (y**2))           = '',1PD15.8)') sum

!     check integration of u(x,y,z) = 0.33*d(y**3)/dy

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(2,j)**3;  v(3) = 3*coor(2,j)**2;  v(2) = 0.d0
            v(4) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(3,j1,i)*v(1)+phi(3,j1+1,i)*v(2)+phi(3,j1+2,i)*v(3)+ &
                       phi(3,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      sum = sum/3.d0
      WRITE (10,'('' Int (0.33*d(y**3)/dy) = '',1PD15.8)') sum

!     check integration of u(x,y,z) = z**2

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(3,j)**2;  v(4) = 2*coor(3,j);  v(3) = 0.d0
            v(2) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(1,j1,i)*v(1)+phi(1,j1+1,i)*v(2)+phi(1,j1+2,i)*v(3)+ &
                       phi(1,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      WRITE (10,'('' Int (z**2)            = '',1PD15.8)') sum

!     check integration of u(x,y,z) = 0.33*d(z**3)/dz

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(3,j)**3;  v(4) = 3*coor(3,j)**2;  v(3) = 0.d0
            v(2) = 0.d0
            j1 = 4*(j-1)+1
            sum = sum+(phi(4,j1,i)*v(1)+phi(4,j1+1,i)*v(2)+phi(4,j1+2,i)*v(3)+ &
                       phi(4,j1+3,i)*v(4))*weight(i)
         ENDDO
      ENDDO
      sum = sum/3d0
      WRITE (10,'('' Int (0.33*d(z**3)/dz) = '',1PD15.8)') sum

!     next integration scheme

      IF (npi.EQ.5) THEN
         npi = 11;  GOTO 10
      ELSEIF (npi.EQ.11) THEN
         npi = 14;  GOTO 10
      ELSEIF (npi.EQ.14) THEN
         npi = 15;  GOTO 10
      ELSEIF (npi.EQ.15) THEN
         npi = 24;  GOTO 10
      ELSEIF (npi.EQ.24) THEN
         npi = 31;  GOTO 10
      ENDIF

      CLOSE (UNIT=10)         
      END


