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

!     We compute numerically the integrals of some functions for which we can 
!     get the exact solution. 
!     From the values of those 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).

!     This does only check the values of the shape functions (stored into phi(1,...)), 
!     not their derivatives stored into phi(2-4,...).

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

      DIMENSION coor(3,4),point(3,14),phi(4,16,14),weight(14),v(4),cu(3,4),cu1(3,4)
      REAL*4 xx

      DATA cu /0.d0,0.d0,0.d0, 1.d0,0.d0,0.d0, 0.d0,1.d0,0.d0, 0.d0,0.d0,1.d0/
      DATA cu1/1.d0,0.d0,0.d0, 1.d0,1.d0,0.d0, 0.d0,1.d0,0.d0, 0.d0,0.d0,1.d0/

      OPEN (UNIT=10,file='resu1.txt',STATUS='UNKNOWN')
      WRITE (10,1)
    1 FORMAT (' Hermite (Délèze) tetrahedron interpolation with 16 parameters'/            &
              ' Check integration of quadratic terms.'/)
      kcas = 1
    2 IF (kcas.EQ.1) THEN
         coor = cu
         WRITE (10,3) kcas
    3    FORMAT (/' Case ',i2,': Unit tetrahedron'/ &
                  ' =========================') 
      ELSEIF (kcas.EQ.2) THEN
         coor = cu1
         WRITE (10,4) kcas
    4    FORMAT (/' Case ',i2,': Unit rotated tetrahedron'/ &
                  ' ================================') 
      ELSEIF (kcas.EQ.3) THEN
         iran = 1.
         DO i = 1,3
            DO j = 1,4
               xx = RAND (iran);  iran = 0
               coor(i,j) = xx*10.d0
            ENDDO
         ENDDO
         WRITE (10,5) kcas
    5    FORMAT (/' Case ',i2,': Random tetrahedron'/ &
                  ' ===========================') 
      ENDIF
      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,'('' 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

!     volume
     
      vol = ABS ( (x2-x1)*((y3-y1)*(z4-z1)-(y4-y1)*(z3-z1)) +                  &
                  (x3-x1)*((y4-y1)*(z2-z1)-(y2-y1)*(z4-z1)) +                  &
                  (x4-x1)*((y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)))
      vol = vol/6.d0 
      WRITE (10,'('' Volume   = '',1PD15.6)') vol

!     centroid

      xm = (coor(1,1)+coor(1,2)+coor(1,3)+coor(1,4))/4.d0
      ym = (coor(2,1)+coor(2,2)+coor(2,3)+coor(2,4))/4.d0
      zm = (coor(3,1)+coor(3,2)+coor(3,3)+coor(3,4))/4.d0

!     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) = 1

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = 1.d0;  v(2) = 0.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
      sume = vol
      WRITE (10,*) 'Int (1)              = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(1,j)-xm;  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
      sume = 0.d0
      WRITE (10,*) 'Int (x-xm)           = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(2,j)-ym;  v(2) = 0.d0;  v(3) = 1.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
      sume = 0.d0
      WRITE (10,*) 'Int (y-ym)           = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = coor(3,j)-zm;  v(2) = 0.d0;  v(3) = 0.d0;  v(4) = 1.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
      sume = 0.d0
      WRITE (10,*) 'Int (z-zm)           = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(1,j)-xm)**2;  v(2) = 2*(coor(1,j)-xm);  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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(1,i)-xm)**2
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((x-xm)**2)      = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(2,j)-ym)**2;  v(2) = 0.d0; v(3) = 2*(coor(2,j)-ym);
            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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(2,i)-ym)**2
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((y-ym)**2)      = ',sum
      WRITE (10,*) 'exact                = ',sume

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

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(3,j)-zm)**2;  v(2) = 0.d0; v(3) = 0.d0
            v(4) = 2*(coor(3,j)-zm);
            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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(3,i)-zm)**2
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((z-zm)**2)      = ',sum
      WRITE (10,*) 'exact                = ',sume
   
!     check integration of u(x,y,z) = (x-xm)*(y-ym)

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(1,j)-xm)*(coor(2,j)-ym);  v(2) = (coor(2,j)-ym)
            v(3) = (coor(1,j)-xm);  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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(1,i)-xm)*(coor(2,i)-ym)
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((x-xm)*(y-ym))  = ',sum
      WRITE (10,*) 'exact                = ',sume

!     check integration of u(x,y,z) = (x-xm)*(z-zm)

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(1,j)-xm)*(coor(3,j)-zm);  v(2) = (coor(3,j)-zm)
            v(3) = 0.d0;  v(4) = (coor(1,j)-xm)
            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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(1,i)-xm)*(coor(3,i)-zm)
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((x-xm)*(z-zm))  = ',sum
      WRITE (10,*) 'exact                = ',sume

!     check integration of u(x,y,z) = (y-ym)*(z-zm)

      sum = 0.d0
      DO i = 1,npi
         DO j = 1,4
            v(1) = (coor(2,j)-ym)*(coor(3,j)-zm);  v(2) = 0.d0
            v(3) = (coor(3,j)-zm);  v(4) = (coor(2,j)-ym)
            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
      sume = 0.d0
      DO i = 1,4
         sume = sume+(coor(2,i)-ym)*(coor(3,i)-zm)
      ENDDO
      sume = sume*vol/20.d0
      WRITE (10,*) 'Int ((y-xm)*(z-ym))  = ',sum
      WRITE (10,*) 'exact                = ',sume

!     next integration scheme

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

!     next tetra

      IF (kcas.EQ.1) THEN
         kcas = 2;  GOTO 2
      ELSEIF (kcas.EQ.2) THEN
         kcas = 3;  GOTO 2
      ENDIF

      CLOSE (UNIT=10)         
      END


