      SUBROUTINE FElib_Tetra0_Shape (coor,point,phi,weight,nnode,npi)

!     +------------------------------------------------------------------------+
!     | Compute the coordinates and values of the shape functions and their    |
!     | first derivatives in a 4, 10, 16 or 20 nodes tetrahedron with linear   |
!     | or quadratic or cubic interpolation, respectively.                     |
!     +------------------------------------------------------------------------+
!     | coor       : coordinates of the 4 main nodes.                          |
!     | point      : coordinates of the integration points in the actual       |
!     |              element.                                                  |
!     | phi(1,i,k) : values of the i-th shape function at the k-th point.      |
!     | phi(2,i,k) : values of the x-derivatives.                              |
!     | phi(3,i,j) : values of the y-derivatives.                              |
!     | phi(4,i,j) : values of the z-derivatives.                              |
!     | weight     : numerical integration weight.                             |
!     | nnode      :  4 for linear interpolation;                              |
!     |              10 for quadratic interpolation;                           |
!     |              16 for cubic interpolation;                               |
!     |              20 for cubic interpolation with bubbles.                  |
!     | npi        : nb. of integration points on the reference tetrahedron    |
!     |              (may be 1, 5, 10).                                        |
!     +------------------------------------------------------------------------+
!     | In case nnode = 4 (linear interpolation):                              |
!     | f1 = t1;  f2 = t2;  f3 = t3;  f4 = t4;                                 |
!     | where ti are the volume coordinates.                                   | 
!     +------------------------------------------------------------------------+
!     | In case nnode = 10 (quadratic interpolation), the mid-side nodes are:  |
!     | node 5 on the edge 1-2;  node 6 on the edge 1-3                        |
!     | node 7 on the edge 1-4;  node 8 on the edge 2-3                        |
!     | node 9 on the edge 3-4;  node 10 on the edge 2-4                       |
!     | With this numbering, the shape functions in volume coordinates are:    |
!     | f1 = (2t1-1)t1;  f2 = (2t2-1)t2;  f3 = (2t3-1)t3;  f4 = (2t4-1)t4;     |
!     | f5 = 4t1t2;  f6 = 4t1t3;  f7 = 4t1t4;  f8 = 4t2t3;  f9 = 4t3t4         |
!     | f10 = 4t2t4.                                                           |
!     +------------------------------------------------------------------------+
!     | In case nnode = 16 or 20 (cubic interpolation), the side nodes are:    |
!     | node 5,6   on the edge 1-2;  node 7,8   on the edge 1-3                |
!     | node 9,10  on the edge 1-4;  node 11,12 on the edge 2-3                |
!     | node 13,14 on the edge 3-4;  node 15,16 on the edge 2-4                |
!     | (node 5 closer to 1, 6 closer to 2; etc.)                              |
!     | f5  = 4.5*t1*t2*(3t1-1); f6  = 4.5*t1*t2*(3t2-1);                      |
!     | f7  = 4.5*t1*t3*(3t1-1); f8  = 4.5*t1*t3*(3t3-1);                      |
!     | f9  = 4.5*t1*t4*(3t1-1); f10 = 4.5*t1*t4*(3t4-1);                      |
!     | f11 = 4.5*t2*t3*(3t2-1); f12 = 4.5*t2*t3*(3t3-1);                      |
!     | f13 = 4.5*t3*t4*(3t3-1); f14 = 4.5*t3*t4*(3t4-1);                      |
!     | f15 = 4.5*t2*t4*(3t2-1); f16 = 4.5*t2*t4*(3t4-1).                      |
!     +------------------------------------------------------------------------+
!     | In case nnode = 20 (cubic interpolation with bubble functions),        |
!     | the mid-face nodes are:                                                |
!     | node 17    at the middle of face 1,2,3;                                |
!     | node 18    at the middle of face 1,2,4;                                |
!     | node 19    at the middle of face 1,3,4;                                |
!     | node 20    at the middle of face 2,3,4;                                |
!     | f17 = 27*t1*t2*t3;  f18 = 27*t1*t2*t4;  f19 = 27*t1*t3*t4;             |
!     | f20 = 27*t2*t3*t4.                                                     |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION coor(3,4),point(3,npi),phi(4,nnode,npi),weight(npi)

!     Local variables:

      DIMENSION a(3,3),dt1(3),dt2(3),dt3(3),dt4(4),tt(4)

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

      INCLUDE 'Tetra_Integration_5.inc'
      INCLUDE 'Tetra_Integration_10.inc'
      INCLUDE 'Tetra_Integration_14.inc'
      INCLUDE 'Tetra_Integration_24.inc'

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

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

      a(1,1) = coor(1,2)-coor(1,1);  a(2,1) = coor(2,2)-coor(2,1)
      a(3,1) = coor(3,2)-coor(3,1)
      a(1,2) = coor(1,3)-coor(1,1);  a(2,2) = coor(2,3)-coor(2,1)
      a(3,2) = coor(3,3)-coor(3,1)
      a(1,3) = coor(1,4)-coor(1,1);  a(2,3) = coor(2,4)-coor(2,1)
      a(3,3) = coor(3,4)-coor(3,1)

!     volume (*6) of the tetrahedron

      x = coor(1,1);  y = coor(2,1);  z = coor(3,1)
      i2 = 2;  i3 = 3;  i4 = 4
      CALL Volume_Coordinate ()
      vol = det

!     coordinates of the integration point in the current element frame

      DO k = 1,npi
         DO i = 1,3
            point(i,k) = coor(i,1)
            DO j = 1,3
               IF (npi.EQ.1) THEN
                  point(i,k) = point(i,k)+a(i,j)*0.25d0
               ELSEIF (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.14) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref14(j,k)
               ELSEIF (npi.EQ.24) THEN
                  point(i,k) = point(i,k)+a(i,j)*ref24(j,k)
               ENDIF
            ENDDO
         ENDDO
      ENDDO

!     numerical integration weights

      DO i = 1,npi
         IF (npi.EQ.1) THEN
            weight(1) = vol/6.d0
         ELSEIF (npi.EQ.5) THEN
            weight(i) = w5(i)*vol/6.d0
         ELSEIF (npi.EQ.10) THEN
            weight(i) = w10(i)*vol/6.d0
         ELSEIF (npi.EQ.14) THEN
            weight(i) = w14(i)*vol/6.d0
         ELSEIF (npi.EQ.24) THEN
            weight(i) = w24(i)*vol/6.d0
         ENDIF
      ENDDO

!     compute the functions phi(1-nnode) and their derivatives

      DO k = 1,npi

!        coordinates of the integration point

         x = point(1,k);  y = point(2,k);  z = point(3,k)

!        volume coordinates ti (i=1-4) of the point (x,y,z)
!        and their derivatives dti(1-3) with respect to x,y,z

         i2 = 2;  i3 = 3;  i4 = 4
         CALL Volume_Coordinate ()
         t1 = det/vol;  dt1(1) = detx/vol;  dt1(2) = dety/vol;  dt1(3) = detz/vol
         i2 = 4;  i3 = 3;  i4 = 1
         CALL Volume_Coordinate ()
         t2 = det/vol;  dt2(1) = detx/vol;  dt2(2) = dety/vol;  dt2(3) = detz/vol
         i2 = 4;  i3 = 1;  i4 = 2
         CALL Volume_Coordinate ()
         t3 = det/vol;  dt3(1) = detx/vol;  dt3(2) = dety/vol;  dt3(3) = detz/vol
         i2 = 2;  i3 = 1;  i4 = 3
         CALL Volume_Coordinate ()
         t4 = det/vol;  dt4(1) = detx/vol;  dt4(2) = dety/vol;  dt4(3) = detz/vol
         tt(1) = t1;  tt(2) = t2;  tt(3) = t3;  tt(4) = t4

!        values of the shape functions at (x,y,z)

         IF (nnode.EQ.4) THEN
            phi(1,1,k) = t1
            phi(1,2,k) = t2
            phi(1,3,k) = t3
            phi(1,4,k) = t4
         ELSE
            phi(1,1,k) = (2.d0*t1-1.d0)*t1
            phi(1,2,k) = (2.d0*t2-1.d0)*t2
            phi(1,3,k) = (2.d0*t3-1.d0)*t3
            phi(1,4,k) = (2.d0*t4-1.d0)*t4
            IF (nnode.EQ.10) THEN
               phi(1,5,k)  = 4.d0*t1*t2
               phi(1,6,k)  = 4.d0*t1*t3
               phi(1,7,k)  = 4.d0*t1*t4
               phi(1,8,k)  = 4.d0*t2*t3
               phi(1,9,k)  = 4.d0*t3*t4
               phi(1,10,k) = 4.d0*t2*t4
            ELSE
               phi(1,5,k)  = 4.5d0*t1*t2*(3d0*t1-1.d0) 
               phi(1,6,k)  = 4.5d0*t1*t2*(3d0*t2-1.d0) 
               phi(1,7,k)  = 4.5d0*t1*t3*(3d0*t1-1.d0) 
               phi(1,8,k)  = 4.5d0*t1*t3*(3d0*t3-1.d0) 
               phi(1,9,k)  = 4.5d0*t1*t4*(3d0*t1-1.d0) 
               phi(1,10,k) = 4.5d0*t1*t4*(3d0*t4-1.d0) 
               phi(1,11,k) = 4.5d0*t2*t3*(3d0*t2-1.d0) 
               phi(1,12,k) = 4.5d0*t2*t3*(3d0*t3-1.d0) 
               phi(1,13,k) = 4.5d0*t3*t4*(3d0*t3-1.d0) 
               phi(1,14,k) = 4.5d0*t3*t4*(3d0*t4-1.d0) 
               phi(1,15,k) = 4.5d0*t2*t4*(3d0*t2-1.d0) 
               phi(1,16,k) = 4.5d0*t2*t4*(3d0*t4-1.d0)
               IF (nnode.EQ.20) THEN
                  phi(1,17,k) = 27.d0*t1*t2*t3 
                  phi(1,18,k) = 27.d0*t1*t2*t4 
                  phi(1,19,k) = 27.d0*t1*t3*t4 
                  phi(1,20,k) = 27.d0*t2*t3*t4 
               ENDIF
            ENDIF
         ENDIF

!        values of the first derivatives with respect to x, y, z

         IF (nnode.EQ.4) THEN                    ! linear 4 nodes
            phi(2,1,k) = dt1(1);  phi(3,1,k) = dt1(2);  phi(4,1,k) = dt1(3)
            phi(2,2,k) = dt2(1);  phi(3,2,k) = dt2(2);  phi(4,2,k) = dt2(3)
            phi(2,3,k) = dt3(1);  phi(3,3,k) = dt3(2);  phi(4,3,k) = dt3(3)
            phi(2,4,k) = dt4(1);  phi(3,4,k) = dt4(2);  phi(4,4,k) = dt4(3)
         ELSE
            DO j = 1,3                           ! quadratic / cubic f1-f4
               phi(j+1,1,k) = (4.d0*t1-1.d0)*dt1(j)
               phi(j+1,2,k) = (4.d0*t2-1.d0)*dt2(j)
               phi(j+1,3,k) = (4.d0*t3-1.d0)*dt3(j)
               phi(j+1,4,k) = (4.d0*t4-1.d0)*dt4(j)
            ENDDO
            IF (nnode.EQ.10) THEN                ! quadratic f5-f10
               DO j = 1,3
                  phi(j+1,5,k)  = 4.d0*(t2*dt1(j)+t1*dt2(j))
                  phi(j+1,6,k)  = 4.d0*(t3*dt1(j)+t1*dt3(j))
                  phi(j+1,7,k)  = 4.d0*(t4*dt1(j)+t1*dt4(j))
                  phi(j+1,8,k)  = 4.d0*(t3*dt2(j)+t2*dt3(j))
                  phi(j+1,9,k)  = 4.d0*(t3*dt4(j)+t4*dt3(j))
                  phi(j+1,10,k) = 4.d0*(t2*dt4(j)+t4*dt2(j))
               ENDDO
            ELSE                                 ! cubic f5-f16
               CALL Partial_Derivative (1,2,0,5)
               CALL Partial_Derivative (2,1,0,6)
               CALL Partial_Derivative (1,3,0,7)
               CALL Partial_Derivative (3,1,0,8)
               CALL Partial_Derivative (1,4,0,9)
               CALL Partial_Derivative (4,1,0,10)
               CALL Partial_Derivative (2,3,0,11)
               CALL Partial_Derivative (3,2,0,12)
               CALL Partial_Derivative (3,4,0,13)
               CALL Partial_Derivative (4,3,0,14)
               CALL Partial_Derivative (2,4,0,15)
               CALL Partial_Derivative (4,2,0,16)
               IF (nnode.EQ.20) THEN             ! cubic f17-f20
                  CALL Partial_Derivative (1,2,3,17)
                  CALL Partial_Derivative (1,2,4,18)
                  CALL Partial_Derivative (1,3,4,19)
                  CALL Partial_Derivative (2,3,4,20)
               ENDIF
            ENDIF
         ENDIF
      ENDDO
      CONTAINS
         SUBROUTINE Volume_Coordinate ()

!        +---------------------------------------------------------------------+
!        | Compute the volume of the sub-tetrahedron (x,y,z)-i2-i3-i4 and its  |
!        | derivatives with respect to x,y,z.                                  |
!        +---------------------------------------------------------------------+

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

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

         x21 = coor(1,i2)-x;  y21 = coor(2,i2)-y;  z21 = coor(3,i2)-z
         x31 = coor(1,i3)-x;  y31 = coor(2,i3)-y;  z31 = coor(3,i3)-z
         x41 = coor(1,i4)-x;  y41 = coor(2,i4)-y;  z41 = coor(3,i4)-z
         det = x21*(y31*z41-y41*z31) + y21*(z31*x41-z41*x31) + z21*(x31*y41-x41*y31)

!        derivatives

         detx =  y41*z31-y31*z41 + y21*z41-y41*z21 + y31*z21-y21*z31
         dety = -x41*z31+x31*z41 - x21*z41+x41*z21 - x31*z21+x21*z31
         detz =  x41*y31-x31*y41 + x21*y41-x41*y21 + x31*y21-x21*y31
         END SUBROUTINE Volume_Coordinate

         SUBROUTINE Partial_Derivative (k1,k2,k3,ii)

!        +---------------------------------------------------------------------+
!        |  Derivatives of phi,x phi,y, phi,z from phi,t1 ... phi,t4           |
!        |  phi(ii) = 4.5*tt(k1)*tt(k2)*(3tt(k1)-1) for ii=1-16                |
!        |  phi(ii) = 27*tt(k1)*tt(k2)*tt(k3)       for ii=17-20               |
!        +---------------------------------------------------------------------+

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

         DIMENSION pt(4)

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

         pt = 0.d0
         IF (k3.EQ.0) THEN
            pt(k1) = 4.5d0*tt(k2)*(6.d0*tt(k1)-1.d0)
            pt(k2) = 4.5d0*tt(k1)*(3.d0*tt(k1)-1.d0)
         ELSE
            pt(k1) = 27.d0*tt(k2)*tt(k3)
            pt(k2) = 27.d0*tt(k1)*tt(k3)
            pt(k3) = 27.d0*tt(k1)*tt(k2)
         ENDIF
         DO j = 1,3
            phi(j+1,ii,k) = pt(1)*dt1(j)+pt(2)*dt2(j)+pt(3)*dt3(j)+pt(4)*dt4(j)
         ENDDO
         END SUBROUTINE Partial_Derivative
      END

      SUBROUTINE FElib_Elem101_Stiff (coor,nnode,imate,rmate,npi,stiff,ndof)

!     +------------------------------------------------------------------------+
!     | Create the stiffness matrix of an elastic tetrahedron with Lagrange's  |
!     | linear interpolation, 'nnode' nodes and 3 degrees of freedom per node, |
!     | for a total of 'ndof' dof.                                             |
!     | Shape functions: see FElib_Tetra0_Shape.f95.                           |
!     +------------------------------------------------------------------------+
!     | coor   : nodal coordinates of the four main nodes.                     |
!     | nnode  : nb. of nodes (may be 4, 10, 16, or 20).                       |
!     | imate  : integer material data.                                        |
!     |          (imate(1) = 1 for an elastic material)                        |
!     | rmate  : real material data.                                           |
!     | npi    : nb. of integration points (may be 1, 5, 10).                  |
!     | stiff  : upper half stiffness matrix.                                  |
!     | ndof   : total nb. of degrees of freedom (=3*nnode).                   |
!     +------------------------------------------------------------------------+

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

!     Parameters:

      DIMENSION coor(3,4),imate(*),rmate(*),stiff(ndof,ndof)

!     Local variables:

      REAL*8, DIMENSION(:),     ALLOCATABLE :: weight
      REAL*8, DIMENSION(:,:),   ALLOCATABLE :: point,b
      REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: phi

      DIMENSION d(6,6)

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

!     check input data

      IF (npi.NE.1.AND.npi.NE.5.AND.npi.NE.10.AND.npi.NE.14.AND.npi.NE.24) npi = 5

!     allocate arrays

      ALLOCATE (point(3,npi),phi(4,nnode,npi),weight(npi),b(6,ndof))

!     compute the shape functions 'phi' and weights 'weight' of the tetrahedron

      CALL FElib_Tetra0_Shape (coor,point,phi,weight,nnode,npi)

!     create the strain to stress matrix 'd' : sigma = d * epsilon
!     material is assumed homogeneous within the element;
!     it may, however, be isotropic or not, depending on the data 'imate' & 'rmate'.

      CALL FElib_Mate_Elastic_3D (imate,rmate,d)

!     initialize stiff

      stiff = 0.d0

!     loop on the integration points 

      DO k = 1,npi

!        Displacements to strain matrix b(6,12); variables are numbered as follows:
!        1-3   : u1,v1,w1  node 1
!        4-6   : u2,v2,w2  node 2
!        etc.

         b = 0.d0
         DO i = 1,nnode
            i1 = 3*(i-1)
            b(1,i1+1) = phi(2,i,k)     ! epsi_x = du/dx           : b(1,*)
            b(2,i1+2) = phi(3,i,k)     ! epsi_y = dv/dy           : b(2,*)
            b(3,i1+3) = phi(4,i,k)     ! epsi_z = dw/dz           : b(3,*)
            b(4,i1+1) = phi(3,i,k)     ! gamma_xy = du/dy + dv/dx : b(4,*)
            b(4,i1+2) = phi(2,i,k)
            b(5,i1+2) = phi(4,i,k)     ! gamma_yz = dv/dz + dw/dy : b(5,*)
            b(5,i1+3) = phi(3,i,k)
            b(6,i1+3) = phi(2,i,k)     ! gamma_zx = dw/dx + du/dz : b(6,*)
            b(6,i1+1) = phi(4,i,k)
         ENDDO

!        contribution to the stiffness matrix: stiff = stiff + trans(b)*d*b*weight

         CALL FElib_Matrix_Symm_BtAB_Add (ndof,6,stiff,d,b,weight(k))
      ENDDO

!     return arrays

      DEALLOCATE (point,phi,weight,b)
      END



