!     +------------------------------------------------------------------------+
!     |   TETRAHEDRON FINITE ELEMENT of Hermite type, class C1, degree two,    |
!     |   with 16 degrees of freedom (u, dxu, dyu, dzu at the nodes).          |
!     |   Variant I-C.                                                         |
!     |   Author: Marcel Délèze, Le Chêne 21, 1633 Marsens, Switzerland.       |
!     |   Ph.D. thesis University of Fribourg, February 1978;                  |
!     |   material made available in November 2008.                            |
!     |   (see ElemFinisTetraC1.pdf for details of the derivation)             |
!     +------------------------------------------------------------------------+
!     |   Modification to Fortran 95 style by G. Dupuis, September 2010        |
!     |   - real in double precision                                           |
!     |   - replace GOTO by IF THEN ... ENDIF making the code more readable.   |
!     |   - change real tests like x.EQ.y by ABS(x-y).LT.small.                |
!     |   - remove unused parameters in routines.                              |
!     |   - add the prefix FElib_Tetra1 to the original routine names.         |
!     |   - Fortran keywords in upper case, variables in lower case.           |
!     |   - general reshape in Fortran 95 style.                               |
!     +------------------------------------------------------------------------+
!     |   File Tetra1.f95 of the library FElib.                                |
!     +------------------------------------------------------------------------+

!     List of routines of this file:

!     Element routines:

!     SUBROUTINE FElib_Tetra1_Elem0 (x1,y1,z1, x2,y2,z2, x3,y3,z3, x4,y4,z4)
!     SUBROUTINE FElib_Tetra1_Elem1
!     SUBROUTINE FElib_Tetra1_Elem2
!     SUBROUTINE FElib_Tetra1_Elem3
!     SUBROUTINE FElib_Tetra1_Herm (px,py,pz, v)
!     SUBROUTINE FElib_Tetra1_Dherm (vx,vy,vz, px,py,pz, dv)

!     Internal routines:

!     SUBROUTINE FElib_Tetra1_Du44 (nfct, vx,vy,vz, px,py,pz, du)
!     SUBROUTINE FElib_Tetra1_Dut44 (nfct, tx,ty,tz, xt,yt,zt, dut)
!     SUBROUTINE FElib_Tetra1_Fe2p16 (c, s)
!     SUBROUTINE FElib_Tetra1_Fe2p44 (c, s)
!     SUBROUTINE FElib_Tetra1_U44 (nfct, px,py,pz, u)
!     SUBROUTINE FElib_Tetra1_Ut44 (nfct, xt,yt,zt, ut)

!     REAL*8 FUNCTION FElib_Tetra1_D1nod (vx,vy,vz, m, c)
!     REAL*8 FUNCTION FElib_Tetra1_Dfint (vx,vy,vz, px,py,pz, s)
!     REAL*8 FUNCTION FElib_Tetra1_Dpt1 (j, tx,ty,tz, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Drt10 (tx,ty,tz, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Drt13 (tx,ty,tz, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Dst10 (tx,ty,tz, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Dst13 (tx,ty,tz, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Dwt11 (tx,ty,tz, xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Fint (px,py,pz, s)
!     REAL*8 FUNCTION FElib_Tetra1_Pt1 (j, wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Rt10 (wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Rt13 (wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_St10 (wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_St13 (wt,xt,yt,zt)
!     REAL*8 FUNCTION FElib_Tetra1_Wt11 (xt,yt,zt)

!     Routines common to Tetra1 and Tetra2

!     SUBROUTINE FElib_Tetra_Decomp (nx,nn,a,scales,ips,d)
!     SUBROUTINE FElib_Tetra_Solve (nx,nn,a,b,x,ips)
!     SUBROUTINE FElib_Tetra_Shape (coor,npi,point,ivers,phi,weight)

      SUBROUTINE FElib_Tetra1_Elem0 (x1,y1,z1, x2,y2,z2, x3,y3,z3, x4,y4,z4)

!     +------------------------------------------------------------------------+
!     | Initialization of the tetrahedron parameters.                          |
!     | input   :  (x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4) coordinates  |
!     |            of the vertices.                                            |
!     | output  :  Tetras common block containing:                             |
!     |             ll     matrix relating the current tetra to the unit one,  |
!     |                    factorized according to the Gauss method.           |
!     |             iperm  new numbering of the rows of ll following pivoting. |
!     |             mx(i), my(i), mz(i), normal vertor to the i-th face.       |
!     |             nx(i,j), ny(i,j), nz(i,j), vectors in the i-th faces,      |
!     |                    normal to the tetra edges (i=1,4).                  |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Local variables:

      DIMENSION aux(3)

!     ==========================================================================

      x(1) = x1;  y(1) = y1;  z(1) = z1
      x(2) = x2;  y(2) = y2;  z(2) = z2
      x(3) = x3;  y(3) = y3;  z(3) = z3
      x(4) = x4;  y(4) = y4;  z(4) = z4

!     matrix ll

      ll(1,1) = x(2)-x(1);  ll(2,1) = y(2)-y(1);  ll(3,1) = z(2)-z(1)
      ll(1,2) = x(3)-x(1);  ll(2,2) = y(3)-y(1);  ll(3,2) = z(3)-z(1)
      ll(1,3) = x(4)-x(1);  ll(2,3) = y(4)-y(1);  ll(3,3) = z(4)-z(1)

!     factorization of ll

      CALL FElib_Tetra_Decomp (3,3,ll,aux,iperm,det)

!     vol = abs(det)/6. is the volume of the thetrahedron;
!     if vol = 0, the subroutine Decomp will display an error message.

      DO i = 1,4
         j = MOD(i,4)+1;  k = MOD(j,4)+1;  l = MOD(k,4)+1

!        vectors (mx(i), my(i), mz(i)).

         mx(i) = (y(k)-y(j))*(z(l)-z(j)) - (z(k)-z(j))*(y(l)-y(j))
         my(i) = (z(k)-z(j))*(x(l)-x(j)) - (x(k)-x(j))*(z(l)-z(j))
         mz(i) = (x(k)-x(j))*(y(l)-y(j)) - (y(k)-y(j))*(x(l)-x(j))
         s = SQRT (mx(i)**2 + my(i)**2 + mz(i)**2)
         mx(i) = mx(i)/s;  my(i) = my(i)/s;  mz(i) = mz(i)/s

!        vectors (nx(i,j), ny(i,j), nz(i,j)),  j = 1,2,3.

         nx(i,1) = my(i)*(z(l)-z(k)) - mz(i)*(y(l)-y(k))
         ny(i,1) = mz(i)*(x(l)-x(k)) - mx(i)*(z(l)-z(k))
         nz(i,1) = mx(i)*(y(l)-y(k)) - my(i)*(x(l)-x(k))
         s = SQRT (nx(i,1)**2+ny(i,1)**2+nz(i,1)**2)
         nx(i,1) = nx(i,1)/s;  ny(i,1) = ny(i,1)/s;  nz(i,1) = nz(i,1)/s

         nx(i,2) = my(i)*(z(j)-z(l)) - mz(i)*(y(j)-y(l))
         ny(i,2) = mz(i)*(x(j)-x(l)) - mx(i)*(z(j)-z(l))
         nz(i,2) = mx(i)*(y(j)-y(l)) - my(i)*(x(j)-x(l))
         s = SQRT (nx(i,2)**2+ny(i,2)**2+nz(i,2)**2)
         nx(i,2) = nx(i,2)/s;  ny(i,2) = ny(i,2)/s;  nz(i,2) = nz(i,2)/s

         nx(i,3) = my(i)*(z(k)-z(j)) - mz(i)*(y(k)-y(j))
         ny(i,3) = mz(i)*(x(k)-x(j)) - mx(i)*(z(k)-z(j))
         nz(i,3) = mx(i)*(y(k)-y(j)) - my(i)*(x(k)-x(j))
         s = SQRT (nx(i,3)**2+ny(i,3)**2+nz(i,3)**2)
         nx(i,3) = nx(i,3)/s;  ny(i,3) = ny(i,3)/s;  nz(i,3) = nz(i,3)/s
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Elem1

!     +------------------------------------------------------------------------+
!     | Interpolation on the tetrahedron.                                      |
!     | Initialization of the diagonal blovks of the matrix q.                 |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     ==========================================================================

      DO i = 1,4
         j = MOD(i,4)+1;  k = MOD(j,4)+1;  l = MOD(k,4)+1

!        initialization of d1, d2, d3 et d4.

         d(i,1,1) = x(j)-x(i);  d(i,1,2) = y(j)-y(i);  d(i,1,3) = z(j)-z(i)
         d(i,2,1) = x(k)-x(i);  d(i,2,2) = y(k)-y(i);  d(i,2,3) = z(k)-z(i)
         d(i,3,1) = x(l)-x(i);  d(i,3,2) = y(l)-y(i);  d(i,3,3) = z(l)-z(i)

!        initialization of d5.

         ind = 3*i-2
         d5(ind  ) = nx(i,1)*(2.*x(j)-x(k)-x(l)) +                             &
                     ny(i,1)*(2.*y(j)-y(k)-y(l)) +                             &
                     nz(i,1)*(2.*z(j)-z(k)-z(l))
         d5(ind+1) = nx(i,2)*(2.*x(k)-x(l)-x(j)) +                             &
                     ny(i,2)*(2.*y(k)-y(l)-y(j)) +                             &
                     nz(i,2)*(2.*z(k)-z(l)-z(j))
         d5(ind+2) = nx(i,3)*(2.*x(l)-x(j)-x(k)) +                             &
                     ny(i,3)*(2.*y(l)-y(j)-y(k)) +                             &
                     nz(i,3)*(2.*z(l)-z(j)-z(k))

!        initialization of d6.

         d6(i) = mx(i)*(3.*x(i)-x(j)-x(k)-x(l)) +                              &
                 my(i)*(3.*y(i)-y(j)-y(k)-y(l)) +                              &
                 mz(i)*(3.*z(i)-z(j)-z(k)-z(l))

!    	 initialisation de d7.

         ind = 3*i-2
         d7(ind  ) = d6(i);  d7(ind+1) = d6(i);  d7(ind+2) = d6(i)
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Elem2

!     +------------------------------------------------------------------------+
!     | Interpolation on the tetrahedron.                                      |
!     | Initialization of the off-diagonal blocks of the matrix q.             |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Local variables:

      DIMENSION du(44)

!     ==========================================================================

      DO i = 1,4
         j = MOD(i,4)+1;  k = MOD(j,4)+1;  l = MOD(k,4)+1

!        initialization of q1.

         px = 0.5*(x(k)+x(l));  py = 0.5*(y(k)+y(l));  pz = 0.5*(z(k)+z(l))
         CALL FElib_Tetra1_Du44 (16,nx(i,1),ny(i,1),nz(i,1),px,py,pz,du)
         ind = 3*i-2
         DO m = 1,16
            q1(ind,m) = du(m)
         ENDDO
         px = 0.5*(x(j)+x(l));  py = 0.5*(y(j)+y(l));  pz = 0.5*(z(j)+z(l))
         CALL FElib_Tetra1_Du44 (16,nx(i,2),ny(i,2),nz(i,2),px,py,pz,du)
         ind = ind+1
         DO m = 1,16
            q1(ind,m) = du(m)
         ENDDO
         px = 0.5*(x(k)+x(j));  py = 0.5*(y(k)+y(j));  pz = 0.5*(z(k)+z(j))
         CALL FElib_Tetra1_Du44 (16,nx(i,3),ny(i,3),nz(i,3),px,py,pz,du)
         ind = ind+1
         DO m = 1,16
            q1(ind,m) = du(m)
         ENDDO

!        initialization of q2.

         px = (x(j)+x(k)+x(l))/3.;  py = (y(j)+y(k)+y(l))/3.
         pz = (z(j)+z(k)+z(l))/3.
         CALL FElib_Tetra1_Du44 (28,mx(i),my(i),mz(i),px,py,pz,du)
         DO m = 1,28
            q2(i,m) = du(m)
         ENDDO

!        initialization of q3.

         px = (4.*x(j)+x(k)+x(l))/6.;  py = (4.*y(j)+y(k)+y(l))/6.
         pz = (4.*z(j)+z(k)+z(l))/6.
         CALL FElib_Tetra1_Du44 (32,mx(i),my(i),mz(i),px,py,pz,du)
         ind = 3*i-2
         DO m = 1,32
            q3(ind,m) = du(m)
         ENDDO
         px = (x(j)+4.*x(k)+x(l))/6.;  py = (y(j)+4.*y(k)+y(l))/6.
         pz = (z(j)+4.*z(k)+z(l))/6.
         CALL FElib_Tetra1_Du44 (32,mx(i),my(i),mz(i),px,py,pz,du)
         ind = ind+1
         DO m = 1,32
            q3(ind,m) = du(m)
         ENDDO
         px = (x(j)+x(k)+4.*x(l))/6.;  py = (y(j)+y(k)+4.*y(l))/6.
         pz = (z(j)+z(k)+4.*z(l))/6.
         CALL FElib_Tetra1_Du44 (32,mx(i),my(i),mz(i),px,py,pz,du)
         ind = ind+1
         DO m = 1,32
            q3(ind,m) = du(m)
         ENDDO
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Elem3

!     +------------------------------------------------------------------------+
!     | Interpolation on the tetrahedron.                                      |
!     | Initialization of the matrix t.                                        |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Local variables:

      DIMENSION c(44), s(44)

!     ==========================================================================

      c = 0.D0
      DO i = 1,16

!        compute the coordinates t(1,i), ..., t(44,i) of the element of the 
!        Hermite's base with respect to the enhanced interpolation base.

         c(i) = 1.D0
         CALL FElib_Tetra1_Fe2p16 (c,s)
         DO j = 1,44
            t(j,i) = s(j)
         ENDDO
         c(i) = 0.D0
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Herm (px,py,pz, v)

!     +------------------------------------------------------------------------+
!     | Compute the values of the Hermite functions at a specified point.      |
!     +------------------------------------------------------------------------+
!     | input   :  Tetras common block initialied by Elem0, ..., Elem3.        | 
!     |            px, py, pz coordinates of the point.                        |
!     | output  :  v(1), ..., v(16) values of the 16 functions.                |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Parameters:

      DIMENSION v(*)

!     Local variables:

      DIMENSION u(44)

!     ==========================================================================

      CALL FElib_Tetra1_U44 (44,px,py,pz, u)
      DO i = 1,16
         si = 0.D0
         DO j = 1,44
            si = si + t(j,i)*u(j)
         ENDDO
         v(i) = si
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Dherm (vx,vy,vz, px,py,pz, dv)

!     +------------------------------------------------------------------------+
!     | Compute the values of the derivatives of the Hermite functions at a    |
!     | specified point.                                                       |
!     +------------------------------------------------------------------------+
!     | input   :  Tetras common block initialied by Elem0, ..., Elem3.        | 
!     |            vx,vy,vz   vector defining the direction of the derivative. |
!     |            px,py,pz   coordinates of the point.                        |
!     | output  :  dv(1), ..., dv(16) values of the derivatives of the 16      |
!     |            functions at p in the direction v.                          |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Parameters:

      DIMENSION dv(*)

!     Local variables:

      DIMENSION du(44)

!     ==========================================================================

      CALL FElib_Tetra1_Du44 (44,vx,vy,vz,px,py,pz,du)
      DO i = 1,16
         si = 0.D0
         DO j = 1,44
            si = si + t(j,i)*du(j)
         ENDDO
         dv(i) = si
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Du44 (nfct, vx,vy,vz, px,py,pz, du)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | input   :  Tetras common block initialied by Elem0, ..., Elem3.        | 
!     |            nfct      can take the values 16, 28, 32 or 44.             |
!     |            vx,vy,vz  vector defining the direction of the derivatives. |
!     |            px,py,pz  a point of the tetrahedron.                       |
!     | output  :  du(1), ..., du(nfct) derivatives of the nfct first          |
!     |                                 functions in the direction (vx,vy,vz)  |
!     |                                 at (px, py, pz).                       |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)
!     Parameters:

      DIMENSION du(*)

!     Local variables:

      DIMENSION gen(3), ref(3)

      DATA nerr/0/, eps/-1.D-03/

!     ==========================================================================

      gen(1) = vx;  gen(2) = vy;  gen(3) = vz
      CALL FElib_Tetra_Solve (3,3,ll,gen,ref,iperm)
      tx = ref(1);  ty = ref(2);  tz = ref(3)
      gen(1) = px-x(1);  gen(2) = py-y(1);  gen(3) = pz-z(1)
      CALL FElib_Tetra_Solve (3,3,ll,gen,ref,iperm)
      xt = ref(1);  yt = ref(2);  zt = ref(3)
      IF (xt.GT.eps .AND. yt.GT.eps .AND. zt.GT.eps .AND.                      &
          1.-xt-yt-zt.GT.eps) GOTO 10

!     the point (px,py,pz) is outside the tetrahedron.

      IF (nerr.LT.10) THEN 
         WRITE (unit = *, fmt = *)                                &
         ' *** error message called from FElib_Tetra1_Du44 :',    &
         ' the point is outside the element.'
         WRITE (unit = *, fmt = *) nfct
         WRITE (unit = *, fmt = *) xt,yt,zt
         WRITE (unit = *, fmt = *) px,py,pz
         WRITE (unit = *, fmt = *) vx,vy,vz

      ENDIF
      nerr = nerr+1

!     derivatives of the base functions.

   10 CALL FElib_Tetra1_Dut44 (nfct,tx,ty,tz,xt,yt,zt,du)
      END

      SUBROUTINE FElib_Tetra1_Dut44 (nfct, tx,ty,tz, xt,yt,zt, dut)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | Derivatives of ut(1), ..., ut(nfct) in the direction (tx,ty,tz) at the |
!     | point (xt,yt,zt).                                                      |
!     | nfct can take the values 16, 28, 32 ou 44.                             |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION dut(*)

!     ==========================================================================

      n = nfct
      t1 = tx;  t2 = ty;  t3 = tz;  t0 = -t1-t2-t3
      x = xt;  y = yt;  z = zt;  w = 1.-x-y-z
      DO i = 1,4

!        block 1.   functions nb. 1 to 16.

         ind = 4*i-3
         dut(ind  ) = FElib_Tetra1_Dpt1 (0,t1,t2,t3,w,x,y,z)
         dut(ind+1) = FElib_Tetra1_Dpt1 (1,t1,t2,t3,w,x,y,z)
         dut(ind+2) = FElib_Tetra1_Dpt1 (1,t2,t3,t1,w,y,z,x)
         dut(ind+3) = FElib_Tetra1_Dpt1 (1,t3,t1,t2,w,z,x,y)
         IF (n.LE.16) GOTO 10

!        block 2.   functions nb. 17 to 28.

         ind = 3*i+14
         dut(ind  ) = FElib_Tetra1_Dwt11 (t1,t2,t3,x,y,z)
         dut(ind+1) = FElib_Tetra1_Dwt11 (t2,t3,t1,y,z,x)
         dut(ind+2) = FElib_Tetra1_Dwt11 (t3,t1,t2,z,x,y)
         IF (n.LE.28) GOTO 10

!        block 3.   functions nb. 29 a 32.

         dut(i+28 ) = FElib_Tetra1_Drt10 (t1,t2,t3,w,x,y,z)
         IF (n.LE.32) GOTO 10

!        block 4.   functions nb. 33 to 44.

         ind = 3*i+30
         dut(ind  ) = FElib_Tetra1_Drt13 (t2,t3,t1,w,y,z,x)
         dut(ind+1) = FElib_Tetra1_Drt13 (t3,t1,t2,w,z,x,y)
         dut(ind+2) = FElib_Tetra1_Drt13 (t1,t2,t3,w,x,y,z)

   10    ch = w;  w = x;  x = y;  y = z;  z = ch
         ch = t0;  t0 = t1;  t1 = t2;  t2 = t3;  t3 = ch
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_Fe2p16 (c, s)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | input   : Tetras common block initialied by Elem0, ..., Elem3.         | 
!     |           c(1), ..., c(16) parameters of the element.                  |
!     | output  : s(1), ..., s(44) coordinates of the interpolation functions  |
!     |                            in the interpolation base of the enhanced   |
!     |                            element.                                    |
!     | c(17), ..., c(44) are used as working space.                           |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Parameters:

      DIMENSION c(*), s(*)

!     ==========================================================================

      DO i = 1,4
         j = MOD(i,4)+1;  k = MOD(j,4)+1;  l = MOD(k,4)+1

!        compute the parameters c(17) a c(28).

         ind = 3*i+14
         c(ind  ) = 0.5*( FElib_Tetra1_D1nod (nx(i,1),ny(i,1),nz(i,1),k,c) +   &
                          FElib_Tetra1_D1nod (nx(i,1),ny(i,1),nz(i,1),l,c) )
         c(ind+1) = 0.5*( FElib_Tetra1_D1nod (nx(i,2),ny(i,2),nz(i,2),l,c) +   &
                          FElib_Tetra1_D1nod (nx(i,2),ny(i,2),nz(i,2),j,c) )
         c(ind+2) = 0.5*( FElib_Tetra1_D1nod (nx(i,3),ny(i,3),nz(i,3),j,c) +   &
                          FElib_Tetra1_D1nod (nx(i,3),ny(i,3),nz(i,3),k,c) )

!        compute the parameters c(29) to c(32).

         dj = FElib_Tetra1_D1nod (mx(i),my(i),mz(i),j,c)
         dk = FElib_Tetra1_D1nod (mx(i),my(i),mz(i),k,c)
         dl = FElib_Tetra1_D1nod (mx(i),my(i),mz(i),l,c)
         c(i+28) = (dj+dk+dl)/3.

!        compute the parameters c(33) to c(44).

         ind = 3*i+30
         c(ind  ) = (4.*dj+dk+dl)/6.
         c(ind+1) = (dj+4.*dk+dl)/6.
         c(ind+2) = (dj+dk+4.*dl)/6.
      ENDDO

!     compute the coordinates s(1), ... , s(44).

      CALL FElib_Tetra1_Fe2p44 (c,s)
      END

      SUBROUTINE FElib_Tetra1_Fe2p44 (c, s)

!     +------------------------------------------------------------------------+
!     | Interpolation of the enhanced element.                                 |
!     | input   : Tetras common block initialied by Elem0, ..., Elem3.         | 
!     |           c(1), ..., c(44) element paraeters.                          |
!     | output  : s(1), ..., s(44) coordinates of the interpolation functions  |
!     |                            in the interpolation base of the enhanced   |
!     |                            element.                                    |
!     | method  : solve the linear system q*s = c.                             |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Parameters:

      DIMENSION c(*), s(*)

!     ==========================================================================

!     compute s(1), ..., s(16).

      DO i = 1,4
         k = 4*i-3
         s(k)   = c(k)
         s(k+1) = d(i,1,1)*c(k+1)+d(i,1,2)*c(k+2)+d(i,1,3)*c(k+3)
         s(k+2) = d(i,2,1)*c(k+1)+d(i,2,2)*c(k+2)+d(i,2,3)*c(k+3)
         s(k+3) = d(i,3,1)*c(k+1)+d(i,3,2)*c(k+2)+d(i,3,3)*c(k+3)
      ENDDO

!     compute s(17), ..., s(28).

      DO i = 1,12
         si = 0.D0
         DO j = 1,16
            si = si+q1(i,j)*s(j)
         ENDDO
         k = i+16
         s(k) = d5(i)*(c(k)-si)
      ENDDO

!     compute s(29), ..., s(32).

      DO i = 1,4
         si = 0.D0
         DO j = 1,28
            si = si+q2(i,j)*s(j)
         ENDDO
         k = 28+i
         s(k) = d6(i)*(c(k)-si)
      ENDDO

!     compute s(33), ..., s(44).

      DO i = 1,12
         si = 0.D0
         DO j = 1,32
            si = si+q3(i,j)*s(j)
         ENDDO
         k = i+32
         s(k) = d7(i)*(c(k)-si)
      ENDDO
      END

      SUBROUTINE FElib_Tetra1_U44 (nfct, px,py,pz, u)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | input   :  Tetras common block initialied by Elem0, ..., Elem3.        | 
!     |            nfct  can thake the values 16, 28, 32 or 44.                |
!     |            px, py, pz a point of the tetrahedron.                      |
!     | output  :  u(1), ..., u(nfct) values of the nfct first base functions  |
!     |                               at (px, py, pz).                         |
!     +------------------------------------------------------------------------+

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

      REAL*8 ll, mx,my,mz, nx,ny,nz
      COMMON /Common_Tetras/ x(4),y(4),z(4),                                   &
             nx(4,3),ny(4,3),nz(4,3), mx(4),my(4),mz(4),                       &
             d(4,3,3),d5(12),d6(4),d7(12), q1(12,16),q2(4,28),q3(12,32),       &
             t(44,16), ll(3,3),iperm(3)

!     Parameters:

      DIMENSION u(*)

!     Local variables:

      DIMENSION gen(3), ref(3)

      DATA nerr/0/, eps/-1.D-03/

!     ==========================================================================

      gen(1) = px-x(1);  gen(2) = py-y(1);  gen(3) = pz-z(1)
      CALL FElib_Tetra_Solve (3,3,ll,gen,ref,iperm)
      xt = ref(1);  yt = ref(2);  zt = ref(3)
      IF (xt.GT.eps .AND. yt.GT.eps .AND. zt.GT.eps .AND. 1.-xt-yt-zt.GT.eps) GOTO 10

!     the (px,py,pz) is outside the tetrahedron.

      IF (nerr.LT.10) WRITE (unit = *, fmt = *)                                &
                      ' *** error message called from FElib_Tetra1_U44 :',     &
                      ' the point is outside the element.'
      nerr = nerr+1

!     values of the base functions.

   10 CALL FElib_Tetra1_Ut44 (nfct,xt,yt,zt,u)
      END

      SUBROUTINE FElib_Tetra1_Ut44 (nfct, xt,yt,zt, ut)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Values of ut(1), ... , ut(nfct) at the point (xt,yt,zt) of the         |
!     | reference tetrahedron.                                                 |
!     | nfct can take the values 16, 28, 32 or 44.                             |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION ut(*)

!     ==========================================================================

      n = nfct;  x = xt;  y = yt;  z = zt;  w = 1.-x-y-z
      DO i = 1,4

!        block 1.    functions nb. 1 to 16.

         ind = 4*i - 3
         ut(ind  ) = FElib_Tetra1_Pt1 (0,w,x,y,z)
         ut(ind+1) = FElib_Tetra1_Pt1 (1,w,x,y,z)
         ut(ind+2) = FElib_Tetra1_Pt1 (1,w,y,z,x)
         ut(ind+3) = FElib_Tetra1_Pt1 (1,w,z,x,y)
         IF (n.LE.16) GOTO 10

!        block 2.    functions nb. 17 to 28.

         ind = 3*i+14
         ut(ind  ) = FElib_Tetra1_Wt11 (x,y,z)
         ut(ind+1) = FElib_Tetra1_Wt11 (y,z,x)
         ut(ind+2) = FElib_Tetra1_Wt11 (z,x,y)
         IF (n.LE.28) GOTO 10

!        block 3.    functions nb. 29 to 32.

         ut(i+28 ) = FElib_Tetra1_Rt10 (w,x,y,z)
         IF (n.LE.32) GOTO 10

!        block 4.    functions nb. 33 to 44.

         ind = 3*i+30
         ut(ind  ) = FElib_Tetra1_Rt13 (w,y,z,x)
         ut(ind+1) = FElib_Tetra1_Rt13 (w,z,x,y)
         ut(ind+2) = FElib_Tetra1_Rt13 (w,x,y,z)

   10    ch = w;  w = x;  x = y;  y = z;  z = ch
      ENDDO
      END

      REAL*8 FUNCTION FElib_Tetra1_D1nod (vx,vy,vz, m, c)

!     +------------------------------------------------------------------------+
!     | input  :  m index of a vertex 1, 2, 3 or 4.                            |
!     |           c(4*m-2), c(4*m-1) and c(4*m) partial derivatives with       |
!     |           respect to x,y,z at vertex m.                                |
!     |           vx, vy, vz direction vector.                                 |
!     | output :  d1nod  derivative in the direction v, at the vertex m.       |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION c(*)

!     ==========================================================================

      ind = 4*m-2
      FElib_Tetra1_D1nod = vx*c(ind) + vy*c(ind+1) + vz*c(ind+2)
      END

      REAL*8 FUNCTION FElib_Tetra1_Dfint (vx,vy,vz, px,py,pz, s)

!     +------------------------------------------------------------------------+
!     | input   :  s(1), ..., s(44) coordinates of an interpolation function   |
!     |            of the enhanced element.                                    |
!     |            vx, vy, vz direction vector.                                |
!     |            px, py, pz) point of the tetrahedron.                       |
!     | output  :  dfint derivative of the interpolation function of           |
!     |            coordinates s in the direction v, at p.                     |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION s(*)

!     Local variables:

      DIMENSION du(44)

!     ==========================================================================

      CALL FElib_Tetra1_Du44 (44,vx,vy,vz,px,py,pz,du)
      dv = 0.D0
      DO i = 1,44
         dv = dv+s(i)*du(i)
      ENDDO
      FElib_Tetra1_Dfint = dv
      END

      REAL*8 FUNCTION FElib_Tetra1_Dpt1 (j, tx,ty,tz, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | Derivative of the polynomial of degree three pt(1,j) in the direction  |
!     | (tx,ty,tz) at point (xt,yt,zt).                                        |
!     | j can take the value 0 ou 1;  wt = 1.-xt-yt-zt.                        |
!     +------------------------------------------------------------------------+

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

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      t1 = tx;  t2 = ty;  t3 = tz;  t0 = -t1-t2-t3
      IF (j.EQ.0) THEN
         p = 2.-w - 2.*(x*x+y*y+z*z+x*y+y*z+z*x)
         dp = -t0 - 2.*( t1*(2.*x+y+z)+t2*(2.*y+z+x)+t3*(2.*z+x+y) )
      ELSE
         p = x*(w+0.5*(y+z))
         dp  = t0*x + t1*(w+0.5*(y+z)) + (t2+t3)*0.5*x
      ENDIF
      FElib_Tetra1_Dpt1 = t0*p + w*dp
      END

      REAL*8 FUNCTION FElib_Tetra1_Drt10 (tx,ty,tz, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | Derivative of the piecewise rational function rt(1,0) in the direction |
!     | at point (xt,yt,zt);  wt = 1.-xt-yt-zt.                                |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      IF (x.GT.eps .AND. y.GT.eps .AND. z.GT.eps) GOTO 10

!     in the faces 2, 3 et 4, we have:

      FElib_Tetra1_Drt10 = 0.D0
      RETURN

   10 t1 = tx;  t2 = ty;  t3 = tz
      IF (x.GE.y .AND. z.GE.y) THEN                             ! piece 2.
         ch = x;  x = y;  y = ch;  ch = t1;  t1 = t2; t2 = ch
      ELSEIF (x.GE.z .AND. y.GE.z) THEN                         ! piece 3.
         ch = x;  x = z;  z = ch;  ch = t1;  t1 = t3;  t3 = ch
      ENDIF

!     piece 1.

      IF (w.LE.eps) THEN 

!        in the intersection of piece 1 and the face 1, we have:

         FElib_Tetra1_Drt10 =  (-t1-t2-t3)*3.*x*x
      ELSE

!        in the intersection of piece 1 and inside the tetrahedron, we have: 

         FElib_Tetra1_Drt10 = FElib_Tetra1_Dst10 (t1,t2,t3,w,x,y,z)
      ENDIF
      END

      REAL*8 FUNCTION FElib_Tetra1_Drt13 (tx,ty,tz, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the enhanced element.                         |
!     | Derivative of the piecewise rational function rt(1,3) in the direction |
!     | (tx,ty,tz) at point (xt,yt,zt);  wt = 1.-xt-yt-zt.                     |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      IF (x.GE.z .AND. y.GE.z) THEN                             ! piece 3.
         FElib_Tetra1_Drt13 = 0.D0
         RETURN
      ENDIF
      IF (x.LE.eps .OR. y.LE.eps .OR. z.LE.eps) THEN

!        dans les faces 2, 3 et 4, on a

         FElib_Tetra1_Drt13 = 0.D0
         RETURN
      ENDIF
      t1 = tx;  t2 = ty;  t3 = tz
      IF (x.GE.y .AND. z.GE.y) THEN                             ! piece 2.
         ch = x;  x = y;  y = ch;  ch = t1;  t1 = t2;  t2 = ch
      ENDIF

!     piece 1.

      IF (w.LE.eps) THEN

!        in the intersection of piece 1 and face 1, we have:

         FElib_Tetra1_Drt13 = 4.*(-t1-t2-t3)*x*(z-x)
         RETURN
      ENDIF

!     in the intersection of piece 1 and inside the tetrahedron, we have:

      FElib_Tetra1_Drt13 = FElib_Tetra1_Dst13 (t1,t2,t3,w,x,y,z)
      END

      REAL*8 FUNCTION FElib_Tetra1_Dst10 (tx,ty,tz, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Derivative of the rational function, restricted to piece 1, in the     |
!     | direction (tx,ty,tz) at point (xt,yt,zt);  wt = 1.-xt-yt-zt.           |
!     | It is assumed that the point is not on one of the five edges of        |
!     | piece 1, belonging to the border of the tetrahedron.                   |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      ymx = y-x;  zmx = z-x
      t1 = tx;  t2 = ty;  t3 = tz
      t0 = -t1-t2-t3
      IF (ymx.GT.eps .OR. zmx.GT.eps) GOTO 10

!     on the line x = y = z, derivative is zero.

      v1 = 0.D0
      dv1 = 0.D0
      GOTO 20

!     in the rest of piece 1, value and derivative of a are stored into v1 and dv1.

   10 p   = x*w
      dp  = t1*w + t0*x
      v1  = -p*( 3./(1.-3.*x) + 2./(y+zmx) )
      dv1 = -3.*(dp - p*( -3.*t1/(1.-3.*x) ))/(1.-3.*x)                        &
            -2.*(dp - p*( (t2+t3-t1)/(y+zmx) ))/(y+zmx)

!     value and derivative of b are stored into v2 and dv2.

      p   = 6.*x-2.*w-3.+v1
      v2  = -ymx*zmx*p
      dv2 = -((t2-t1)*zmx+ymx*(t3-t1))*p - ymx*zmx*(6.*t1-2.*t0+dv1)

!     value and derivative of c are stored into v1 and dv1.

      p   = zmx*v2
      dp  = (t3-t1)*v2 + zmx*dv2
      q   = (ymx+zmx)*y*(w+ymx)
      v1  = p/q
      dv1 = (dp - p*( (t2+t3-2.*t1)/(ymx+zmx)+t2/y+(t0+t2-t1)/(w+ymx) ))/q
      p   = ymx*v2
      dp  = (t2-t1)*v2 + ymx*dv2
      q   = (ymx+zmx)*z*(w+zmx)
      v1  = v1 + p/q
      dv1 = dv1 + (dp - p*( (t2+t3-2.*t1)/(ymx+zmx)+t3/z + (t0+t3-t1)/(w+zmx) ))/q

!     value and derivative of d are stored into v2 and dv2.

   20 p   = w*( 2.*(ymx+zmx)+v1 )
      dp  = t0*( 2.*(ymx+zmx)+v1)+ w*( 2.*(t2+t3-2.*t1)+dv1 )
      q   = (1.-3.*x)*(y+zmx)
      v2  = p/q
      dv2 = (dp - p*( -3.*t1/(1.-3.*x)+(t2+t3-t1)/(y+zmx) ))/q

!     derivative of rt(1,0).

      FElib_Tetra1_Dst10 = x*( (t1*2.*w+x*t0)*(3.+v2) + x*w*dv2 )
      END

      REAL*8 FUNCTION FElib_Tetra1_Dst13 (tx,ty,tz, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Derivative of the rational function rt(1,3), restricted to piece 1,    |
!     | in the direction (tx,ty,tz) at point (xt,yt,zt);  wt = 1.-xt-yt-zt.    |
!     | It is assumed that the point is not on one of the five edges of        |
!     | piece 1, belonging to the border of the tetrahedron.                   |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt;  ymx = y-x;  zmx = z-x
      IF (ymx.GT.eps .OR. zmx.GT.eps) GOTO 10

!     on the line x = y = z, we have:

      FElib_Tetra1_Dst13 = 0.D0
      RETURN

!     in the rest of piece 1, we have:

   10 t1 = tx;  t2 = ty;  t3 = tz;  t0 = -t1-t2-t3

!     value and derivative of e.

      p  = x*zmx
      dp = t1*zmx + x*(t3-t1)
      q  = 1.-3.*x
      e  = p/q
      de = (dp - e*(-3.*t1))/q
      q  = w+y
      e  = -4.*x+2.*z + e - p/q
      de = -4.*t1+2.*t3 + de - ( dp - p*(t0+t2)/q )/q

!     value and derivative of f.

      p  = w*ymx*e
      dp = t0*ymx*e + w*(t2-t1)*e + w*ymx*de
      q  = y*(ymx+zmx)*(w+ymx)*(w+y)*(1.-3.*x)
      f  = p/q
      df = (dp - p*( t2/y+(t2+t3-2.*t1)/(ymx+zmx)+(t0+t2-t1)/(w+ymx)           &
                 +(t0+t2)/(w+y)-3.*t1/(1.-3.*x) ))/q
      q  = (w+x)*(w+zmx)
      f  = f + 2./q
      df = df - 2.*( (t0+t1)/(w+x)+(t0+t3-t1)/(w+zmx) )/q

!     derivative of dst13.

      FElib_Tetra1_Dst13 = 2.*x*zmx*( 2.*(t1*zmx+x*(t3-t1))*w*f + x*zmx*(t0*f+w*df) )
      END

      REAL*8 FUNCTION FElib_Tetra1_Dwt11 (tx,ty,tz, xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Derivative of the piecewise polynomial wt(1,1) in the direction        |
!     | (tx,ty,tz) at point (xt,yt,zt).                                        |
!     +------------------------------------------------------------------------+

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

!     ==========================================================================

      x = xt;  y = yt;  z = zt
      IF (x.GE.y .AND. z.GE.y) THEN                   ! piece 2.
         FElib_Tetra1_Dwt11 = y*(ty*(2.*z-y) + tz*y)
      ELSEIF (x.GE.z .AND. y.GE.z) THEN               ! piece 3.
         FElib_Tetra1_Dwt11 = z*(ty*z + tz*(2.*y-z))
      ELSE                                            ! piece 1.
         FElib_Tetra1_Dwt11 = 2.*tx*(x*(x-y-z)+y*z) + x*(ty*(2.*z-x)+tz*(2.*y-x))
      ENDIF
      END

      REAL*8 FUNCTION FElib_Tetra1_Fint (px,py,pz, s)

!     +------------------------------------------------------------------------+
!     | input  : s(1), ..., s(44) coordinates of the interlation function in   |
!     |          the base of the enhanced element.                             |
!     |          (px, py, pz) a point of the tetrahedron.                      |
!     | output : fint value of the interpolation function of coordinates s     |
!     |          at point (px, py, pz).                                        |
!     +------------------------------------------------------------------------+

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

!     Paremeters:

      DIMENSION s(*)

!     Local variables:

      DIMENSION u(44)

!     ==========================================================================

      CALL FElib_Tetra1_U44 (44,px,py,pz,u)
      v = 0.D0
      DO i = 1,44
         v = v+s(i)*u(i)
      ENDDO
      FElib_Tetra1_Fint = v
      END

      REAL*8 FUNCTION FElib_Tetra1_Pt1 (j, wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the polynomial of degree three pt(1,j) at point (xt,yt,zt).   |
!     | j can take the value 0 ou 1;  wt = 1.-xt-yt-zt.                        |
!     +------------------------------------------------------------------------+

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

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      IF (j .EQ. 0) THEN
         FElib_Tetra1_Pt1 = w*(2.-w-2.*(x*x+y*y+z*z+x*y+y*z+z*x))
      ELSE
         FElib_Tetra1_Pt1 = x*w*(w+0.5*(y+z))
      ENDIF
      END

      REAL*8 FUNCTION FElib_Tetra1_Rt10 (wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the piecewise rational function rt(1,0) at the point          |
!     | (xt,yt,zt) of the reference element;  wt = 1.-xt-yt-zt.                |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      IF (x.GT.eps .AND. y.GT.eps .AND. z.GT.eps .AND. w.GT.eps ) GOTO 10

!     on the border of the tetrahedron, we have:

      FElib_Tetra1_Rt10 = 0.D0
      RETURN

!     inside the tetrahedron, we have:

   10 IF (x.GT.y .AND. z.GE.y) THEN         ! piece 2.
         ch = x;  x = y;  y = ch
      ELSEIF (x.GE.z .AND. y.GE.z) THEN     ! piece 3.
         ch = x;  x = z;  z = ch
      ENDIF

!     piece 1.

      FElib_Tetra1_Rt10 = FElib_Tetra1_St10 (w,x,y,z)
      END

      REAL*8 FUNCTION FElib_Tetra1_Rt13 (wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the piecewise rational function rt(1,3) at the point          |
!     | (xt,yt,zt) of the reference element;  wt = 1.-xt-yt-zt.                |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      IF (x.GE.z .AND. y.GE.z) THEN              ! piece 3.
         FElib_Tetra1_Rt13 = 0.D0
         RETURN
      ENDIF
      IF (x.GT.eps .AND. y.GT.eps .AND. z.GT.eps .AND. w.GT.eps ) GOTO 10

!     on the border of the tetrahedron, we have:

      FElib_Tetra1_Rt13 = 0.D0
      RETURN

!     inside the tetrahedron, we have:

   10 IF (x.GE.y .AND. z.GE.y) THEN              ! piece 2.
         ch = x;  x = y;  y = ch
      ENDIF

!     piece 1.

      FElib_Tetra1_Rt13 = FElib_Tetra1_St13 (w,x,y,z)
      END

      REAL*8 FUNCTION FElib_Tetra1_St10 (wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the rational function rt(1,0), restricted to piece 1, at the  |
!     | point (xt,yt,zt);  wt = 1.-xt-yt-zt.                                   |
!     | It is assumed that the point is not on one of the five edges of        |
!     | piece 1, belonging to the border of the tetrahedron.                   |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      ymx = y-x;  zmx = z-x
      IF (ymx.LE.eps .AND. zmx.LE.eps) THEN      ! on the line x = y = z
         FElib_Tetra1_St10 = 3.*x*x*w
      ELSE                                       ! in the rest of piece 1
         var = -x*w*( 3./(1.-3.*x) + 2./(y+zmx) )
         var = -ymx*zmx*( 6.*x-2.*w-3.+var )
         var = var*( zmx/(y*(w+ymx))+ymx/(z*(w+zmx)) )/(ymx+zmx)
         var = w*( 2.*(ymx+zmx)+var )/((1.-3.*x)*(y+zmx))
         FElib_Tetra1_St10 = x*x*w*(3.+var)
      ENDIF
      END

      REAL*8 FUNCTION FElib_Tetra1_St13 (wt,xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the rational function rt(1,3), restricted to piece 1, at the  |
!     | point (xt,yt,zt);  wt = 1.-xt-yt-zt.                                   |
!     | It is assumed that the point is not on one of the five edges of        |
!     | piece 1, belonging to the border of the tetrahedron.                   |
!     +------------------------------------------------------------------------+

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

      DATA eps/0.5D-06/

!     ==========================================================================

      w = wt;  x = xt;  y = yt;  z = zt
      ymx = y-x;  zmx = z-x
      IF (ymx.GT.eps .OR. zmx.GT.eps) GOTO 10

!     on the line x = y = z, we have:

      FElib_Tetra1_St13 = 0.D0
      RETURN

!     in the rest of piece 1, we have:

   10 var = 1./(1.-3.*x) - 1./(w+y)
      var = -4.*x + 2.*z + x*zmx*var
      var = w*ymx*var /( y*(ymx+zmx)*(w+ymx)*(w+y)*(1.-3.*x) )
      var = var + 2./( (w+x)*(w+zmx) )
      FElib_Tetra1_St13 = 2.*x*x*zmx*zmx*w*var
      END

      REAL*8 FUNCTION FElib_Tetra1_Wt11 (xt,yt,zt)

!     +------------------------------------------------------------------------+
!     | Base of interpolation of the reference element.                        |
!     | Value of the piecewise polynomial wt(1,1) at point (xt,yt,zt).         |
!     +------------------------------------------------------------------------+

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

!     ==========================================================================

      x = xt;  y = yt;  z = zt
      IF (x.GE.y .AND. z.GE.y) THEN
         FElib_Tetra1_Wt11 = y*y*( z-y/3. )                   ! piece 2.
      ELSEIF (x.GE.z .AND. y.GE.z) THEN
         FElib_Tetra1_Wt11 = z*z*( y-z/3. )                   ! piece 3.
      ELSE
         FElib_Tetra1_Wt11 = x*( 2.*x*x/3.-x*y-x*z+2.*y*z )   ! piece 1.
      ENDIF
      END

      SUBROUTINE FElib_Tetra_Decomp (nx,nn,a,scales,ips,d)

!     +------------------------------------------------------------------------+ 
!     | This routine does an LU decomposition of a square matrix.              |
!     | It uses row scaling and row pivoting.                                  |
!     | The parameters of the calling sequence are:                            |
!     | nn      the order of the matrix to be decomposed.                      |
!     | nx      the number of lines of the matrix a in the calling program.    |
!     | a       the matrix to be decomposed.                                   |
!     | d       the determinant of the coefficient matrix.                     |
!     | kerr    error print control and return error flag                      |
!     |         if kerr is set negative by the user no message is printed      |
!     |         if set to zero, error message will be printed.                 |
!     |         on return, following the code:                                 |
!     |         1 matrix with zero row.                                        |
!     |         2 singular matrix.                                             |
!     |         3 solution obtained.                                           |
!     | ips    vector containing permutated row numbers of a.                  |
!     +------------------------------------------------------------------------+
 
      IMPLICIT REAL*8 (a-h, o-z), INTEGER*4 (i-n)

      DIMENSION a(nx,*), scales(*), ips(*)

      DATA small /1.e-20/

!     ==========================================================================

      n = nn;  kerr = 0
 
!     find the solution in the case n = 1

      IF (n.EQ.1) THEN
         d = a(1,1)
         ips(1) = 1
         IF (ABS(d).LT.small) THEN
            IF (kerr.EQ.0) CALL FElib_Tetra_Sing (1)
            kerr = 1
            RETURN
         ENDIF
         kerr = 3
         RETURN
      ENDIF
 
!     initialize

      odd = 1.
 
!     find the scaling factor for each row and store in scales

      DO i = 1,n
         ips(i) = i
         rmax = 0.D0
         DO j = 1,n
            rmax = MAX (rmax,ABS (a(i,j)))
         ENDDO
         IF (ABS(rmax).LT.small) GOTO 10
         scales(i) = 1./rmax
      ENDDO
 
!     main loop. Gaussian elimination is carried out. Loop is gone
!     through once for each element on the main diagonal. k is the
!     diagonal element number. Last diagonal element is an exception.

      nm1 = n-1
      DO k = 1,nm1
         big = 0.D0
         DO i = k,n
            ip = ips(i)
            size = ABS (a(ip,k))*scales(ip)
            IF (big.GT.size) CYCLE
            big = size;  ipiv = i
         ENDDO
         IF (ABS(big).LT.small) GOTO 20
         IF (ipiv.NE.k) THEN
            j = ips(k);  ips(k) = ips(ipiv);  ips(ipiv) = j;  odd = -odd
         ENDIF
         kp = ips(k);  pivot = a(kp,k);  kp1 = k+1
         DO i = kp1,n
            ip = ips(i);  em = a(ip,k)/pivot;  a(ip,k) = em
            IF (ABS(em).LT.small) CYCLE
            DO j = kp1,n
               a(ip,j) = a(ip,j)-em*a(kp,j)
            ENDDO
         ENDDO
      ENDDO
      kp = ips(n)
      IF (ABS(a(kp,n)).LT.small) GOTO 20
 
!     calculate the determinant

      d = 1.
      DO i = 1,n
         ip = ips(i);  d = d*a(ip,i)
      ENDDO
      d = d*odd
      kerr = 3
      RETURN
 
!     error exist

   10 IF (kerr.EQ.0) CALL FElib_Tetra_Sing (1)
      d = 0.D0
      kerr = 1
      RETURN

   20 IF (kerr.EQ.0) CALL FElib_Tetra_Sing (2)
      d = 0.D0
      kerr = 2
      RETURN
      END
 
      SUBROUTINE FElib_Tetra_Sing (iwhy)
 
!     +------------------------------------------------------------------------+ 
!     | Display an error message.                                              |
!     +------------------------------------------------------------------------+
 
      IF (iwhy.EQ.1) THEN
         WRITE (unit = *, fmt = *) ' *** matrix with zero row in Decomp'
      ELSEIF (iwhy.EQ.2) THEN
         WRITE (unit = *, fmt = *) ' *** singular matrix in Decomp.'
         WRITE (unit = *, fmt = *) '     zero divide in solve.'
      ENDIF
      END
 
      SUBROUTINE FElib_Tetra_Solve (nx,nn,a,b,x,ips)

!     +------------------------------------------------------------------------+ 
!     | This routine does forward elimination and back substitution to solve   |
!     | the equation lu*x = b.                                                 |
!     | The parameters of the calling sequence are:                            |
!     | nn      the order of lu.                                               |
!     | nx      the number of lines of the matrix lu in the calling program.   |
!     | a       the decomposed matrix, l and u are packed into one matrix by   |
!     |         Decomp.                                                        |
!     | b       the right hand side of the equation.                           |
!     | ips     vector containing the permutated row numbers.                  |
!     +------------------------------------------------------------------------+ 

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

      DIMENSION a(nx,*), b(*), x(*), ips(*)

!     ==========================================================================

      n = nn
 
!     find the solution in the case n = 1

      IF (n.EQ.1) THEN
         x(1) = b(1)/a(1,1)
         RETURN
      ENDIF
 
!     forward elimination using l and storing the result into x

      np1 = n+1;  ip = ips(1);  x(1) = b(ip)
      DO i = 2,n
         ip = ips(i);  im1 = i-1;  sum = 0.D0
         DO j = 1,im1
            sum = sum+a(ip,j)*x(j)
         ENDDO
         x(i) = b(ip)-sum
      ENDDO
 
!     back substitution using u and x

      ip = ips(n);  x(n) = x(n)/a(ip,n)
      DO iback = 2,n
         i = np1-iback;  ip = ips(i);  ip1 = i+1;  sum = 0.D0
         DO j = ip1,n
            sum = sum+a(ip,j)*x(j)
         ENDDO
         x(i) = (x(i)-sum)/a(ip,i)
      ENDDO
      END

      SUBROUTINE FElib_Tetra_Shape (coor,npi,point,ivers,phi,weight)

!     +-----------------------------------------------------------------------+
!     | Compute the shape functions and their derivatives of an tetrahedron   |
!     | with 16 degrees of freedom and Délèze's interpolation.                |
!     +-----------------------------------------------------------------------+
!     | coor  : nodal coordinates (coor(3,4))                                 |
!     | npi   : nb. of integration points (minimum is 5)                      |
!     | point : coordinates of the integration points (point(3,npi))          |
!     | ivers : type of interpolation (1 or 2)                                |
!     | phi   : shape functions phi(i,j,k) of phi(4,16,npi)                   |
!     |         i = 1 for the function                                        |
!     |             2 x-derivative; 3 y-derivative; 4 z-derivative.           |
!     |         j = function index as follows:                                |
!     |             phi1,dxphi1,dyphi1,dzphi1, ..., phi4,dxphi4,dyphi4,dzphi4 |
!     |         k = index of the integration point                            |
!     | weight: integration weights related to the npi points.                |
!     +-----------------------------------------------------------------------+

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

!     Parameters:
   
      DIMENSION coor(3,4),point(3,*),phi(4,16,*),weight(*)

!     Local variables:

      DIMENSION a(3,3),u(16),orig(3),w5(5),ref5(3,5),w10(10),ref10(3,10),      &
                w11(11),ref11(3,11),w14(14),ref14(3,14),w15(15),ref15(3,15),   &
                w24(24),ref24(3,24),w31(31),ref31(3,31) 

!     coordinates of the integration point in the reference element frame
!     and integration weights

!     5 points precision 2

      DATA w5    /-0.8000000000000000, 0.4500000000000000, 0.4500000000000000, &
                   0.4500000000000000, 0.4500000000000000/
      DATA ref5  / 0.2500000000000000, 0.2500000000000000, 0.2500000000000000, &
                   0.5000000000000000, 0.1666666666666667, 0.1666666666666667, &
                   0.1666666666666667, 0.1666666666666667, 0.1666666666666667, &
                   0.1666666666666667, 0.1666666666666667, 0.5000000000000000, &
                   0.1666666666666667, 0.5000000000000000, 0.1666666666666667/

!     10 points precision 3

      DATA w10   / 0.2177650698804054, 0.2177650698804054, 0.2177650698804054, &
                   0.2177650698804054, 0.0214899534130631, 0.0214899534130631, &
                   0.0214899534130631, 0.0214899534130631, 0.0214899534130631, &
                   0.0214899534130631/
      DATA ref10 / 0.5684305841968444, 0.1438564719343852, 0.1438564719343852, &
                   0.1438564719343852, 0.1438564719343852, 0.1438564719343852, &
                   0.1438564719343852, 0.1438564719343852, 0.5684305841968444, &
                   0.1438564719343852, 0.5684305841968444, 0.1438564719343852, &
                   0.0000000000000000, 0.5000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.0000000000000000, 0.5000000000000000/

!     11 points precision 4

      DATA w11   /-0.0789333333333333, 0.0457333333333333, 0.0457333333333333, &
                   0.0457333333333333, 0.0457333333333333, 0.1493333333333333, &
                   0.1493333333333333, 0.1493333333333333, 0.1493333333333333, &
                   0.1493333333333333, 0.1493333333333333/ 
      DATA ref11 / 0.2500000000000000, 0.2500000000000000, 0.2500000000000000, &
                   0.7857142857142857, 0.0714285714285714, 0.0714285714285714, &
                   0.0714285714285714, 0.0714285714285714, 0.0714285714285714, &
                   0.0714285714285714, 0.0714285714285714, 0.7857142857142857, &
                   0.0714285714285714, 0.7857142857142857, 0.0714285714285714, &
                   0.1005964238332008, 0.3994035761667992, 0.3994035761667992, &
                   0.3994035761667992, 0.1005964238332008, 0.3994035761667992, &
                   0.3994035761667992, 0.3994035761667992, 0.1005964238332008, &
                   0.3994035761667992, 0.1005964238332008, 0.1005964238332008, &
                   0.1005964238332008, 0.3994035761667992, 0.1005964238332008, &
                   0.1005964238332008, 0.1005964238332008, 0.3994035761667992/

!     14 points precision 5

      DATA w14   / 0.0190476190476190, 0.0190476190476190, 0.0190476190476190, &
                   0.0190476190476190, 0.0190476190476190, 0.0190476190476190, &
                   0.0885898247429807, 0.0885898247429807, 0.0885898247429807, &
                   0.0885898247429807, 0.1328387466855907, 0.1328387466855907, &
                   0.1328387466855907, 0.1328387466855907/
      DATA ref14 / 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.0000000000000000, 0.5000000000000000, &
                   0.6984197043243866, 0.1005267652252045, 0.1005267652252045, &
                   0.1005267652252045, 0.1005267652252045, 0.1005267652252045, &
                   0.1005267652252045, 0.1005267652252045, 0.6984197043243866, &
                   0.1005267652252045, 0.6984197043243866, 0.1005267652252045, &
                   0.0568813795204234, 0.3143728734931922, 0.3143728734931922, &
                   0.3143728734931922, 0.3143728734931922, 0.3143728734931922, &
                   0.3143728734931922, 0.3143728734931922, 0.0568813795204234, &
                   0.3143728734931922, 0.0568813795204234, 0.3143728734931922/

!     15 points precision 5

      DATA w15   / 0.1817020685825351, 0.0361607142857143, 0.0361607142857143, &
                   0.0361607142857143, 0.0361607142857143, 0.0698714945161738, &
                   0.0698714945161738, 0.0698714945161738, 0.0698714945161738, &
                   0.0656948493683187, 0.0656948493683187, 0.0656948493683187, &
                   0.0656948493683187, 0.0656948493683187, 0.0656948493683187/
      DATA ref15 / 0.2500000000000000, 0.2500000000000000, 0.2500000000000000, &
                   0.0000000000000000, 0.3333333333333333, 0.3333333333333333, &
                   0.3333333333333333, 0.3333333333333333, 0.3333333333333333, &
                   0.3333333333333333, 0.3333333333333333, 0.0000000000000000, &
                   0.3333333333333333, 0.0000000000000000, 0.3333333333333333, &
                   0.7272727272727273, 0.0909090909090909, 0.0909090909090909, &
                   0.0909090909090909, 0.0909090909090909, 0.0909090909090909, &
                   0.0909090909090909, 0.0909090909090909, 0.7272727272727273, &
                   0.0909090909090909, 0.7272727272727273, 0.0909090909090909, &
                   0.4334498464263357, 0.0665501535736643, 0.0665501535736643, &
                   0.0665501535736643, 0.4334498464263357, 0.0665501535736643, &
                   0.0665501535736643, 0.0665501535736643, 0.4334498464263357, &
                   0.0665501535736643, 0.4334498464263357, 0.4334498464263357, &
                   0.4334498464263357, 0.0665501535736643, 0.4334498464263357, &
                   0.4334498464263357, 0.4334498464263357, 0.0665501535736643/

!     24 points precision 6

      DATA w24 /   0.0399227502581679, 0.0399227502581679, 0.0399227502581679, &
                   0.0399227502581679, 0.0100772110553207, 0.0100772110553207, &
                   0.0100772110553207, 0.0100772110553207, 0.0553571815436544, &
                   0.0553571815436544, 0.0553571815436544, 0.0553571815436544, &
                   0.0482142857142857, 0.0482142857142857, 0.0482142857142857, &
                   0.0482142857142857, 0.0482142857142857, 0.0482142857142857, &
                   0.0482142857142857, 0.0482142857142857, 0.0482142857142857, &
                   0.0482142857142857, 0.0482142857142857, 0.0482142857142857/ 
      DATA ref24 / 0.3561913862225449, 0.2146028712591517, 0.2146028712591517, & 
                   0.2146028712591517, 0.2146028712591517, 0.2146028712591517, &
                   0.2146028712591517, 0.2146028712591517, 0.3561913862225449, &
                   0.2146028712591517, 0.3561913862225449, 0.2146028712591517, &
                   0.8779781243961660, 0.0406739585346113, 0.0406739585346113, &
                   0.0406739585346113, 0.0406739585346113, 0.0406739585346113, &
                   0.0406739585346113, 0.0406739585346113, 0.8779781243961660, &
                   0.0406739585346113, 0.8779781243961660, 0.0406739585346113, &
                   0.0329863295731731, 0.3223378901422757, 0.3223378901422757, &
                   0.3223378901422757, 0.3223378901422757, 0.3223378901422757, &
                   0.3223378901422757, 0.3223378901422757, 0.0329863295731731, &
                   0.3223378901422757, 0.0329863295731731, 0.3223378901422757, &
                   0.2696723314583159, 0.0636610018750175, 0.0636610018750175, &
                   0.0636610018750175, 0.2696723314583159, 0.0636610018750175, &
                   0.0636610018750175, 0.0636610018750175, 0.2696723314583159, &
                   0.6030056647916491, 0.0636610018750175, 0.0636610018750175, &
                   0.0636610018750175, 0.6030056647916491, 0.0636610018750175, &
                   0.0636610018750175, 0.0636610018750175, 0.6030056647916491, &
                   0.0636610018750175, 0.2696723314583159, 0.6030056647916491, &
                   0.2696723314583159, 0.6030056647916491, 0.0636610018750175, &
                   0.6030056647916491, 0.0636610018750175, 0.2696723314583159, &
                   0.0636610018750175, 0.6030056647916491, 0.2696723314583159, &
                   0.2696723314583159, 0.0636610018750175, 0.6030056647916491, &
                   0.6030056647916491, 0.2696723314583159, 0.0636610018750175/

!     31 points precision 7

      DATA w31 /   0.1095853407966528, 0.0635996491464850, 0.0635996491464850, &
                   0.0635996491464850, 0.0635996491464850,-0.3751064406859797, &
                  -0.3751064406859797,-0.3751064406859797,-0.3751064406859797, &
                   0.0293485515784412, 0.0293485515784412, 0.0293485515784412, &
                   0.0293485515784412, 0.0058201058201058, 0.0058201058201058, &
                   0.0058201058201058, 0.0058201058201058, 0.0058201058201058, &
                   0.0058201058201058, 0.1653439153439105, 0.1653439153439105, &
                   0.1653439153439105, 0.1653439153439105, 0.1653439153439105, &
                   0.1653439153439105, 0.1653439153439105, 0.1653439153439105, &
                   0.1653439153439105, 0.1653439153439105, 0.1653439153439105, &
                   0.1653439153439105/
      DATA ref31 / 0.2500000000000000, 0.2500000000000000, 0.2500000000000000, &
                   0.7653604230090441, 0.0782131923303186, 0.0782131923303186, &
                   0.0782131923303186, 0.0782131923303186, 0.0782131923303186, &
                   0.0782131923303186, 0.0782131923303186, 0.7653604230090441, &
                   0.0782131923303186, 0.7653604230090441, 0.0782131923303186, &
                   0.6344703500082868, 0.1218432166639044, 0.1218432166639044, &
                   0.1218432166639044, 0.1218432166639044, 0.1218432166639044, &
                   0.1218432166639044, 0.1218432166639044, 0.6344703500082868, &
                   0.1218432166639044, 0.6344703500082868, 0.1218432166639044, &
                   0.0023825066607383, 0.3325391644464206, 0.3325391644464206, &
                   0.3325391644464206, 0.3325391644464206, 0.3325391644464206, &
                   0.3325391644464206, 0.3325391644464206, 0.0023825066607383, &
                   0.3325391644464206, 0.0023825066607383, 0.3325391644464206, &
                   0.0000000000000000, 0.5000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.5000000000000000, &
                   0.5000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.5000000000000000, 0.0000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.5000000000000000, 0.0000000000000000, &
                   0.0000000000000000, 0.0000000000000000, 0.5000000000000000, &
                   0.2000000000000000, 0.1000000000000000, 0.1000000000000000, &
                   0.1000000000000000, 0.2000000000000000, 0.1000000000000000, &
                   0.1000000000000000, 0.1000000000000000, 0.2000000000000000, &
                   0.6000000000000000, 0.1000000000000000, 0.1000000000000000, &
                   0.1000000000000000, 0.6000000000000000, 0.1000000000000000, &
                   0.1000000000000000, 0.1000000000000000, 0.6000000000000000, &
                   0.1000000000000000, 0.2000000000000000, 0.6000000000000000, &
                   0.2000000000000000, 0.6000000000000000, 0.1000000000000000, &
                   0.6000000000000000, 0.1000000000000000, 0.2000000000000000, &
                   0.1000000000000000, 0.6000000000000000, 0.2000000000000000, &
                   0.2000000000000000, 0.1000000000000000, 0.6000000000000000, &
                   0.6000000000000000, 0.2000000000000000, 0.1000000000000000/

!     ==========================================================================

!     default parameters

      IF (ivers.NE.1.AND.ivers.NE.2) ivers = 1
      IF (npi.NE. 5.AND.npi.NE.10.AND.npi.NE.11.AND.npi.NE.14.AND.              &
          npi.NE.15.AND.npi.NE.24.AND.npi.NE.31) npi = 5

!     nodal coordinates

      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)

!     mapping of the tetrahedron to the unit (reference) one

      a(1,1) = x2-x1;  a(2,1) = y2-y1;  a(3,1) = z2-z1
      a(1,2) = x3-x1;  a(2,2) = y3-y1;  a(3,2) = z3-z1
      a(1,3) = x4-x1;  a(2,3) = y4-y1;  a(3,3) = z4-z1
      orig(1) = x1;  orig(2) = y1;  orig(3) = z1 
      rjac = ABS ( a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) +                      &
                   a(1,2)*(a(2,3)*a(3,1)-a(2,1)*a(3,3)) +                      &
                   a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1)) )

!     coordinates of the integration point in the current element frame

      DO k = 1,npi
         DO i = 1,3
            point(i,k) = orig(i)
            DO j = 1,3
               IF (npi.EQ.5) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref5(j,k)
               ELSEIF (npi.EQ.10) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref10(j,k)
               ELSEIF (npi.EQ.11) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref11(j,k)
               ELSEIF (npi.EQ.14) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref14(j,k)
               ELSEIF (npi.EQ.15) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref15(j,k)
               ELSEIF (npi.EQ.24) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref24(j,k)
               ELSEIF (npi.EQ.31) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref31(j,k)
               ENDIF
            ENDDO
         ENDDO
      ENDDO

!     numerical integration weights

      IF (npi.EQ.5) THEN
         DO i = 1,npi
            weight(i) = w5(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.10) THEN
         DO i = 1,npi
            weight(i) = w10(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.11) THEN
         DO i = 1,npi
            weight(i) = w11(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.14) THEN
         DO i = 1,npi
            weight(i) = w14(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.15) THEN
         DO i = 1,npi
            weight(i) = w15(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.24) THEN
         DO i = 1,npi
            weight(i) = w24(i)*rjac/6.d0
         ENDDO
      ELSEIF (npi.EQ.31) THEN
         DO i = 1,npi
            weight(i) = w31(i)*rjac/6.d0
         ENDDO
      ENDIF

!     tetra initialization (see Tetra1 and Tetra2 (not yet))

!!      IF (ivers.EQ.1) THEN
         CALL FElib_Tetra1_Elem0 (x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)
         CALL FElib_Tetra1_Elem1 ()
         CALL FElib_Tetra1_Elem2 ()
         CALL FElib_Tetra1_Elem3 ()
!!      ELSE
!!         CALL FElib_Tetra2_Elem0 (x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)
!!         CALL FElib_Tetra2_Elem1 ()
!!         CALL FElib_Tetra2_Elem2 ()
!!         CALL FElib_Tetra2_Elem3 ()
!!      ENDIF 

!     find the values of the shape functions

!!      IF (ivers.EQ.1) THEN
         DO k = 1,npi
            CALL FElib_Tetra1_Herm (point(1,k),point(2,k),point(3,k), u) 
            DO j = 1,16
               phi(1,j,k) = u(j)
            ENDDO
         ENDDO
!!      ELSE
!!         DO k = 1,npi
!!            CALL FElib_Tetra2_Herm (point(1,k),point(2,k),point(3,k), u) 
!!            DO j = 1,16
!!               phi(1,j,k) = u(j)
!!            ENDDO
!!         ENDDO
!!      ENDIF

!     find the values of the derivatives

!!      IF (ivers.EQ.1) THEN
         DO k = 1,npi
            p1 = point(1,k);  p2 = point(2,k);  p3 = point(3,k);  
            dx = 1.d0;  dy = 0.d0;  dz = 0.d0
            CALL FElib_Tetra1_Dherm (dx,dy,dz,p1,p2,p3,u) 
            DO j = 1,16
               phi(2,j,k) = u(j)
            ENDDO
            dx = 0.d0;  dy = 1.d0;  dz = 0.d0
            CALL FElib_Tetra1_Dherm (dx,dy,dz,p1,p2,p3,u)  
            DO j = 1,16
               phi(3,j,k) = u(j)
            ENDDO
            dx = 0.d0;  dy = 0.d0;  dz = 1.d0
            CALL FElib_Tetra1_Dherm (dx,dy,dz,p1,p2,p3,u)  
            DO j = 1,16
               phi(4,j,k) = u(j)
            ENDDO
         ENDDO
!!      ELSE
!!         DO k = 1,npi
!!            p1 = point(1,k);  p2 = point(2,k);  p3 = point(3,k);  
!!            dx = 1.d0;  dy = 0.d0;  dz = 0.d0
!!            CALL FElib_Tetra2_Dherm (dx,dy,dz,p1,p2,p3,u) 
!!            DO j = 1,16
!!               phi(2,j,k) = u(j)
!!            ENDDO
!!            dx = 0.d0;  dy = 1.d0;  dz = 0.d0
!!            CALL FElib_Tetra2_Dherm (dx,dy,dz,p1,p2,p3,u)  
!!            DO j = 1,16
!!               phi(3,j,k) = u(j)
!!            ENDDO
!!            dx = 0.d0;  dy = 0.d0;  dz = 1.d0
!!            CALL FElib_Tetra2_Dherm (dx,dy,dz,p1,p2,p3,u)  
!!            DO j = 1,16
!!               phi(4,j,k) = u(j)
!!            ENDDO
!!         ENDDO
!!      ENDIF

      END

