v3d10.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE v3d10(coor,matcstet,lmcstet,R_in,d,ci, &
     s11l,s22l,s33l,s12l,s23l,s13l,&
     numnp,nstart,nend,numcstet,numat_vol)

!!****f* Rocfrac/Rocfrac/Source/v3d10.f90
!!
!!  NAME
!!    v3d10
!!
!!  FUNCTION
!!
!!    Computes the internal force for a small deformation, 
!!    linear elastic 10-node tetrahedral.
!!
!!  INPUTS
!!
!!   NumNP    -- Number of nodes
!!   numcstet -- Number of elements
!!   Coor     -- number of coordinates
!!   Matcstet  -- Material id
!!   d -- nodal displacement
!!   ci -- compliance matrix
!! 
!!       stored in the follow format:
!!
!!            [ ci(1) ci(2) ci(4)                   ]
!!            [       ci(3) ci(5)                   ]
!!            [             ci(6)                   ]
!!            [                   ci(7)             ]
!!            [                         ci(8)       ]
!!            [                               ci(9) ]
!!
!!   lmcstet  -- Nodal connectivity
!!   nstart, nend -- element beginning and end loop counter
!!   numat_vol -- number of materials
!!
!!  OUTPUT
!!
!!    Rin -- internal force vector
!!    S11l,S22l,S33l,S12l,S23l,S13l -- Stresses at gauss points
!!
!!****

!----- 4 Guass points 
  IMPLICIT NONE
!-----Global variables
  INTEGER :: numnp          ! number of nodes
  INTEGER :: numat_vol      ! number of volumetric materials
  INTEGER :: numcstet       ! number of LSTets
!--   coordinate array
  REAL*8, DIMENSION(1:3,1:numnp) :: coor
!--   elastic stiffness consts
  REAL*8, DIMENSION(1:9,1:numat_vol) :: ci
!--   internal force
  REAL*8, DIMENSION(1:3*numnp) :: r_in
!--   displacement vector
  REAL*8, DIMENSION(1:3*numnp) :: d
!--   CSTet stress
  REAL*8, DIMENSION(1:4,1:numcstet) :: s11l, s22l, s33l, s12l, s23l, s13l
!--   connectivity table for CSTet  
  INTEGER, DIMENSION(1:10,1:numcstet) :: lmcstet
!--   mat number for CSTet element
  INTEGER, DIMENSION(1:numcstet) :: matcstet
!---- Local variables
!--   node numbers
  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
!--   x, y, and z displacements of nodes  
  REAL*8 :: u1,u2,u3,u4,u5,u6,u7,u8,u9,u10
  REAL*8 :: v1,v2,v3,v4,v5,v6,v7,v8,v9,v10
  REAL*8 :: w1,w2,w3,w4,w5,w6,w7,w8,w9,w10
!--   6*volume and inverse of 6*volume      
  REAL*8 :: vx6, vx6inv
!--   spacial derivatives
  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
  REAL*8 :: b11,b12,b13,b14,b15,b16,b17,b18,b19,b20
  REAL*8 :: b21,b22,b23,b24,b25,b26,b27,b28,b29,b30
!--   strains
  REAL*8 :: e11,e22,e33,e12,e23,e13
!--   coordinate holding variable
  REAL*8 :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
  REAL*8 :: y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
  REAL*8 :: z1,z2,z3,z4,z5,z6,z7,z8,z9,z10
!--   dummy and counters
  INTEGER :: i,j,nstart,nend
  REAL*8 :: aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12
!--   partial internal force
  REAL*8 :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,r17,r18
  REAL*8 :: r19,r20,r21,r22,r23,r24,r25,r26,r27,r28,r29,r30
  REAL*8 :: g1, g2, g3, g4
  REAL*8 :: xn1, xn2, xn3, xn4
!--  Coordinate subtractions
  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!--  
  REAL*8 :: c11, c21, c31
  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8,k1n9,k1n10
  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8,k2n9,k2n10
  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8,k3n9,k3n10
  
  
  DO i = nstart, nend
     j   = matcstet(i)
     
     n1  = lmcstet(1,i)
     n2  = lmcstet(2,i)
     n3  = lmcstet(3,i)
     n4  = lmcstet(4,i)
     n5  = lmcstet(5,i)
     n6  = lmcstet(6,i)
     n7  = lmcstet(7,i)
     n8  = lmcstet(8,i)
     n9  = lmcstet(9,i)
     n10 = lmcstet(10,i)
     
     k3n1  = 3*n1
     k3n2  = 3*n2
     k3n3  = 3*n3
     k3n4  = 3*n4
     k3n5  = 3*n5
     k3n6  = 3*n6
     k3n7  = 3*n7
     k3n8  = 3*n8
     k3n9  = 3*n9
     k3n10 = 3*n10
     
     k2n1  = k3n1  - 1
     k2n2  = k3n2  - 1
     k2n3  = k3n3  - 1
     k2n4  = k3n4  - 1
     k2n5  = k3n5  - 1
     k2n6  = k3n6  - 1
     k2n7  = k3n7  - 1
     k2n8  = k3n8  - 1
     k2n9  = k3n9  - 1
     k2n10 = k3n10 - 1
     
     k1n1  = k3n1  - 2
     k1n2  = k3n2  - 2
     k1n3  = k3n3  - 2
     k1n4  = k3n4  - 2 
     k1n5  = k3n5  - 2
     k1n6  = k3n6  - 2
     k1n7  = k3n7  - 2
     k1n8  = k3n8  - 2 
     k1n9  = k3n9  - 2
     k1n10 = k3n10 - 2
                                ! k#n# dummy variables replaces:
     u1  = d(k1n1)          ! (3*n1 -2)
     u2  = d(k1n2)          ! (3*n2 -2)
     u3  = d(k1n3)          ! (3*n3 -2)
     u4  = d(k1n4)          ! (3*n4 -2)
     u5  = d(k1n5)          ! (3*n5 -2)
     u6  = d(k1n6)          ! (3*n6 -2)
     u7  = d(k1n7)          ! (3*n7 -2)
     u8  = d(k1n8)          ! (3*n8 -2)
     u9  = d(k1n9)          ! (3*n9 -2)
     u10 = d(k1n10)         ! (3*n10-2)         
     v1  = d(k2n1)          ! (3*n1 -1)
     v2  = d(k2n2)          ! (3*n2 -1)
     v3  = d(k2n3)          ! (3*n3 -1)
     v4  = d(k2n4)          ! (3*n4 -1)
     v5  = d(k2n5)          ! (3*n5 -1)
     v6  = d(k2n6)          ! (3*n6 -1)
     v7  = d(k2n7)          ! (3*n7 -1)
     v8  = d(k2n8)          ! (3*n8 -1)
     v9  = d(k2n9)          ! (3*n9 -1)
     v10 = d(k2n10)         ! (3*n10-1)
     w1  = d(k3n1)          ! (3*n1)
     w2  = d(k3n2)          ! (3*n2)
     w3  = d(k3n3)          ! (3*n3)
     w4  = d(k3n4)          ! (3*n4)
     w5  = d(k3n5)          ! (3*n5)
     w6  = d(k3n6)          ! (3*n6)
     w7  = d(k3n7)          ! (3*n7)
     w8  = d(k3n8)          ! (3*n8)
     w9  = d(k3n9)          ! (3*n9)
     w10 = d(k3n10)         ! (3*n10)         
     
     x1 = coor(1,n1)
     x2 = coor(1,n2)
     x3 = coor(1,n3)
     x4 = coor(1,n4)
     y1 = coor(2,n1)
     y2 = coor(2,n2)
     y3 = coor(2,n3)
     y4 = coor(2,n4)
     z1 = coor(3,n1)
     z2 = coor(3,n2)
     z3 = coor(3,n3)
     z4 = coor(3,n4)
     
     x14 = x1 - x4
     x24 = x2 - x4
     x34 = x3 - x4
     y14 = y1 - y4
     y24 = y2 - y4
     y34 = y3 - y4
     z14 = z1 - z4
     z24 = z2 - z4
     z34 = z3 - z4
     
     c11 =    y24*z34 - z24*y34
     c21 = -( x24*z34 - z24*x34 )
     c31 =    x24*y34 - y24*x34

     vx6 = -( x14*c11 + y14*c21 + z14*c31 )

!        Vx6 = x2*y3*z4 - x2*y4*z3 - y2*x3*z4 + y2*x4*z3 + z2*x3*y4 
!     $        - z2*x4*y3 - x1*y3*z4 + x1*y4*z3 + x1*y2*z4 - x1*y2*z3
!     $        - x1*z2*y4 + x1*z2*y3 + y1*x3*z4 - y1*x4*z3 - y1*x2*z4 
!     $        + y1*x2*z3 + y1*z2*x4 - y1*z2*x3 - z1*x3*y4 + z1*x4*y3
!     $        + z1*x2*y4 - z1*x2*y3 - z1*y2*x4 + z1*y2*x3

     vx6inv = 1.d0 / vx6
     
     aux1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3)
     aux2 =  (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3)
     aux3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3)
     aux4 =  (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3)
     aux5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3)
     aux6 =  (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3)
     aux7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2)
     aux8 =  (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2)
     aux9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2)
     aux10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2)
     aux11 =-(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2)
     aux12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2)

!--  Gauss Point 1

     g1 = 0.58541020d0
     g2 = 0.13819660d0
     g3 = 0.13819660d0
     g4 = 0.13819660d0

     xn1 = (4.d0*g1-1.d0)   ! derivative of shape function
     xn2 = (4.d0*g2-1.d0)   ! dN_i/dzeta_i
     xn3 = (4.d0*g3-1.d0)
     xn4 = (4.d0*g4-1.d0)
!        xN5 = 4.d0*g1*g2
!        xN6 = 4.d0*g2*g3
!        xN7 = 4.d0*g3*g1
!        xN8 = 4.d0*g1*g4
!        xN9 = 4.d0*g2*g4
!        xN10= 4.d0*g3*g4

     b1 = aux1*xn1
     b2 = aux2*xn1
     b3 = aux3*xn1
     b4 = aux4*xn2
     b5 = aux5*xn2
     b6 = aux6*xn2
     b7 = aux7*xn3
     b8 = aux8*xn3
     b9 = aux9*xn3
     b10 =  aux10*xn4
     b11 =  aux11*xn4
     b12 =  aux12*xn4

     b13 =  4.d0*(g2*aux1 + g1*aux4)
     b14 =  4.d0*(g2*aux2 + g1*aux5)
     b15 =  4.d0*(g2*aux3 + g1*aux6)
     
     b16 =  4.d0*(g3*aux4 + g2*aux7)
     b17 =  4.d0*(g3*aux5 + g2*aux8)
     b18 =  4.d0*(g3*aux6 + g2*aux9)
     
     b19 =  4.d0*(g1*aux7 + g3*aux1)
     b20 =  4.d0*(g1*aux8 + g3*aux2)
     b21 =  4.d0*(g1*aux9 + g3*aux3)
     
     b22 =  4.d0*(g4*aux1 + g1*aux10)
     b23 =  4.d0*(g4*aux2 + g1*aux11)
     b24 =  4.d0*(g4*aux3 + g1*aux12)
     
     b25 =  4.d0*(g4*aux4 + g2*aux10)
     b26 =  4.d0*(g4*aux5 + g2*aux11)
     b27 =  4.d0*(g4*aux6 + g2*aux12)
     
     b28 =  4.d0*(g4*aux7 + g3*aux10)
     b29 =  4.d0*(g4*aux8 + g3*aux11)
     b30 =  4.d0*(g4*aux9 + g3*aux12)
     
     e11 = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + &
             b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10) * vx6inv
     e22 = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + &
          b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10) * vx6inv
     e33 = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + &
          b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10) * vx6inv
     e12 = (b2*u1 + b1*v1 + b5*u2 + b4*v2  &
          + b8*u3 + b7*v3 + b11*u4 + b10*v4 &
          + b14*u5 + b13*v5 + b17*u6 + b16*v6 &
          + b20*u7 + b19*v7 + b23*u8 + b22*v8 &
          + b26*u9 + b25*v9 + b29*u10 + b28*v10) * vx6inv
     e23 = (b3*v1 + b2*w1 + b6*v2 + b5*w2  &
          + b9*v3 + b8*w3 + b12*v4 + b11*w4 &
          + b15*v5 + b14*w5 + b18*v6 + b17*w6 &
          + b21*v7 + b20*w7 + b24*v8 + b23*w8 &
          + b27*v9 + b26*w9 + b30*v10 + b29*w10) * vx6inv
     e13 = (b3*u1 + b1*w1 + b6*u2 + b4*w2  &
          + b9*u3 + b7*w3 + b12*u4 + b10*w4 &
          + b15*u5 + b13*w5 + b18*u6 + b16*w6 &
          + b21*u7 + b19*w7 + b24*u8 + b22*w8 &
          + b27*u9 + b25*w9 + b30*u10 + b28*w10) * vx6inv 
     
     s11l(1,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
     s22l(1,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
     s33l(1,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
     s12l(1,i) = e12*ci(7,j)
     s23l(1,i) = e23*ci(8,j)
     s13l(1,i) = e13*ci(9,j)
     
     r1 = s11l(1,i)*b1 + s12l(1,i)*b2 + s13l(1,i)*b3
     r2 = s22l(1,i)*b2 + s12l(1,i)*b1 + s23l(1,i)*b3
     r3 = s33l(1,i)*b3 + s23l(1,i)*b2 + s13l(1,i)*b1
     
     r4 = s11l(1,i)*b4 + s12l(1,i)*b5 + s13l(1,i)*b6
     r5 = s22l(1,i)*b5 + s12l(1,i)*b4 + s23l(1,i)*b6
     r6 = s33l(1,i)*b6 + s23l(1,i)*b5 + s13l(1,i)*b4
     
     r7 = s11l(1,i)*b7 + s12l(1,i)*b8 + s13l(1,i)*b9
     r8 = s22l(1,i)*b8 + s12l(1,i)*b7 + s23l(1,i)*b9
     r9 = s33l(1,i)*b9 + s23l(1,i)*b8 + s13l(1,i)*b7
     
     r10 = s11l(1,i)*b10 + s12l(1,i)*b11 + s13l(1,i)*b12
     r11 = s22l(1,i)*b11 + s12l(1,i)*b10 + s23l(1,i)*b12
     r12 = s33l(1,i)*b12 + s23l(1,i)*b11 + s13l(1,i)*b10
     
     r13 = s11l(1,i)*b13 + s12l(1,i)*b14 + s13l(1,i)*b15
     r14 = s22l(1,i)*b14 + s12l(1,i)*b13 + s23l(1,i)*b15
     r15 = s33l(1,i)*b15 + s23l(1,i)*b14 + s13l(1,i)*b13
     
     r16 = s11l(1,i)*b16 + s12l(1,i)*b17 + s13l(1,i)*b18
     r17 = s22l(1,i)*b17 + s12l(1,i)*b16 + s23l(1,i)*b18
     r18 = s33l(1,i)*b18 + s23l(1,i)*b17 + s13l(1,i)*b16
     
     r19 = s11l(1,i)*b19 + s12l(1,i)*b20 + s13l(1,i)*b21
     r20 = s22l(1,i)*b20 + s12l(1,i)*b19 + s23l(1,i)*b21
     r21 = s33l(1,i)*b21 + s23l(1,i)*b20 + s13l(1,i)*b19
     
     r22 = s11l(1,i)*b22 + s12l(1,i)*b23 + s13l(1,i)*b24
     r23 = s22l(1,i)*b23 + s12l(1,i)*b22 + s23l(1,i)*b24
     r24 = s33l(1,i)*b24 + s23l(1,i)*b23 + s13l(1,i)*b22
     
     r25 = s11l(1,i)*b25 + s12l(1,i)*b26 + s13l(1,i)*b27
     r26 = s22l(1,i)*b26 + s12l(1,i)*b25 + s23l(1,i)*b27
     r27 = s33l(1,i)*b27 + s23l(1,i)*b26 + s13l(1,i)*b25
     
     r28 = s11l(1,i)*b28 + s12l(1,i)*b29 + s13l(1,i)*b30
     r29 = s22l(1,i)*b29 + s12l(1,i)*b28 + s23l(1,i)*b30
     r30 = s33l(1,i)*b30 + s23l(1,i)*b29 + s13l(1,i)*b28
!--  Gauss Point 2

     g1 = 0.13819660d0
     g2 = 0.58541020d0
     g3 = 0.13819660d0
     g4 = 0.13819660d0

     xn1 = (4.d0*g1-1.d0)   ! derivative of shape function
     xn2 = (4.d0*g2-1.d0)   ! dN_i/dzeta_i
     xn3 = (4.d0*g3-1.d0)
     xn4 = (4.d0*g4-1.d0)
!        xN5 = 4.d0*g1*g2
!        xN6 = 4.d0*g2*g3
!        xN7 = 4.d0*g3*g1
!        xN8 = 4.d0*g1*g4
!        xN9 = 4.d0*g2*g4
!        xN10= 4.d0*g3*g4
     
     b1 = aux1*xn1
     b2 = aux2*xn1
     b3 = aux3*xn1
     b4 = aux4*xn2
     b5 = aux5*xn2
     b6 = aux6*xn2
     b7 = aux7*xn3
     b8 = aux8*xn3
     b9 = aux9*xn3
     b10 =  aux10*xn4
     b11 =  aux11*xn4
     b12 =  aux12*xn4
     b13 =  4.d0*(g2*aux1 + g1*aux4)
     b14 =  4.d0*(g2*aux2 + g1*aux5)
     b15 =  4.d0*(g2*aux3 + g1*aux6)
     b16 =  4.d0*(g3*aux4 + g2*aux7)
     b17 =  4.d0*(g3*aux5 + g2*aux8)
     b18 =  4.d0*(g3*aux6 + g2*aux9)
     b19 =  4.d0*(g1*aux7 + g3*aux1)
     b20 =  4.d0*(g1*aux8 + g3*aux2)
     b21 =  4.d0*(g1*aux9 + g3*aux3)
     b22 =  4.d0*(g4*aux1 + g1*aux10)
     b23 =  4.d0*(g4*aux2 + g1*aux11)
     b24 =  4.d0*(g4*aux3 + g1*aux12)
     b25 =  4.d0*(g4*aux4 + g2*aux10)
     b26 =  4.d0*(g4*aux5 + g2*aux11)
     b27 =  4.d0*(g4*aux6 + g2*aux12)
     b28 =  4.d0*(g4*aux7 + g3*aux10)
     b29 =  4.d0*(g4*aux8 + g3*aux11)
     b30 =  4.d0*(g4*aux9 + g3*aux12)
     
     e11 = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + &
          b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10) * vx6inv 
     e22 = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + &
          b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10) * vx6inv 
     e33 = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + &
          b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10) * vx6inv
     e12 = (b2*u1 + b1*v1 + b5*u2 + b4*v2  &
          + b8*u3 + b7*v3 + b11*u4 + b10*v4 &
          + b14*u5 + b13*v5 + b17*u6 + b16*v6 &
          + b20*u7 + b19*v7 + b23*u8 + b22*v8 &
          + b26*u9 + b25*v9 + b29*u10 + b28*v10) * vx6inv
     e23 = (b3*v1 + b2*w1 + b6*v2 + b5*w2  &
          + b9*v3 + b8*w3 + b12*v4 + b11*w4 &
          + b15*v5 + b14*w5 + b18*v6 + b17*w6 &
          + b21*v7 + b20*w7 + b24*v8 + b23*w8 &
          + b27*v9 + b26*w9 + b30*v10 + b29*w10) * vx6inv
     e13 = (b3*u1 + b1*w1 + b6*u2 + b4*w2  &
          + b9*u3 + b7*w3 + b12*u4 + b10*w4 &
          + b15*u5 + b13*w5 + b18*u6 + b16*w6 &
          + b21*u7 + b19*w7 + b24*u8 + b22*w8 &
          + b27*u9 + b25*w9 + b30*u10 + b28*w10) * vx6inv

     s11l(2,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
     s22l(2,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
     s33l(2,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
     s12l(2,i) = e12*ci(7,j)
     s23l(2,i) = e23*ci(8,j)
     s13l(2,i) = e13*ci(9,j)
     
     r1 =  r1 +s11l(2,i)*b1 + s12l(2,i)*b2 + s13l(2,i)*b3
     r2 =  r2 +s22l(2,i)*b2 + s12l(2,i)*b1 + s23l(2,i)*b3
     r3 =  r3 +s33l(2,i)*b3 + s23l(2,i)*b2 + s13l(2,i)*b1
     r4 =  r4 +s11l(2,i)*b4 + s12l(2,i)*b5 + s13l(2,i)*b6
     r5 =  r5 +s22l(2,i)*b5 + s12l(2,i)*b4 + s23l(2,i)*b6
     r6 =  r6 +s33l(2,i)*b6 + s23l(2,i)*b5 + s13l(2,i)*b4
     r7 =  r7 +s11l(2,i)*b7 + s12l(2,i)*b8 + s13l(2,i)*b9
     r8 =  r8 +s22l(2,i)*b8 + s12l(2,i)*b7 + s23l(2,i)*b9
     r9 =  r9 +s33l(2,i)*b9 + s23l(2,i)*b8 + s13l(2,i)*b7
     r10 = r10 +s11l(2,i)*b10 + s12l(2,i)*b11 + s13l(2,i)*b12
     r11 = r11 +s22l(2,i)*b11 + s12l(2,i)*b10 + s23l(2,i)*b12
     r12 = r12 +s33l(2,i)*b12 + s23l(2,i)*b11 + s13l(2,i)*b10
     r13 = r13 +s11l(2,i)*b13 + s12l(2,i)*b14 + s13l(2,i)*b15
     r14 = r14 +s22l(2,i)*b14 + s12l(2,i)*b13 + s23l(2,i)*b15
     r15 = r15 +s33l(2,i)*b15 + s23l(2,i)*b14 + s13l(2,i)*b13
     r16 = r16 +s11l(2,i)*b16 + s12l(2,i)*b17 + s13l(2,i)*b18
     r17 = r17 +s22l(2,i)*b17 + s12l(2,i)*b16 + s23l(2,i)*b18
     r18 = r18 +s33l(2,i)*b18 + s23l(2,i)*b17 + s13l(2,i)*b16
     r19 = r19 +s11l(2,i)*b19 + s12l(2,i)*b20 + s13l(2,i)*b21
     r20 = r20 +s22l(2,i)*b20 + s12l(2,i)*b19 + s23l(2,i)*b21
     r21 = r21 +s33l(2,i)*b21 + s23l(2,i)*b20 + s13l(2,i)*b19
     r22 = r22 +s11l(2,i)*b22 + s12l(2,i)*b23 + s13l(2,i)*b24
     r23 = r23 +s22l(2,i)*b23 + s12l(2,i)*b22 + s23l(2,i)*b24
     r24 = r24 +s33l(2,i)*b24 + s23l(2,i)*b23 + s13l(2,i)*b22
     r25 = r25 +s11l(2,i)*b25 + s12l(2,i)*b26 + s13l(2,i)*b27
     r26 = r26 +s22l(2,i)*b26 + s12l(2,i)*b25 + s23l(2,i)*b27
     r27 = r27 +s33l(2,i)*b27 + s23l(2,i)*b26 + s13l(2,i)*b25
     r28 = r28 +s11l(2,i)*b28 + s12l(2,i)*b29 + s13l(2,i)*b30
     r29 = r29 +s22l(2,i)*b29 + s12l(2,i)*b28 + s23l(2,i)*b30
     r30 = r30 +s33l(2,i)*b30 + s23l(2,i)*b29 + s13l(2,i)*b28
!--  Gauss Point 3

     g1 = 0.13819660d0
     g2 = 0.13819660d0
     g3 = 0.58541020d0
     g4 = 0.13819660d0
     
     xn1 = (4.d0*g1-1.d0)   ! derivative of shape function
     xn2 = (4.d0*g2-1.d0)   ! dN_i/dzeta_i
     xn3 = (4.d0*g3-1.d0)
     xn4 = (4.d0*g4-1.d0)
!        xN5 = 4.d0*g1*g2
!        xN6 = 4.d0*g2*g3
!        xN7 = 4.d0*g3*g1
!        xN8 = 4.d0*g1*g4
!        xN9 = 4.d0*g2*g4
!        xN10= 4.d0*g3*g4

     b1 = aux1*xn1
     b2 = aux2*xn1
     b3 = aux3*xn1
     b4 = aux4*xn2
     b5 = aux5*xn2
     b6 = aux6*xn2
     b7 = aux7*xn3
     b8 = aux8*xn3
     b9 = aux9*xn3
     b10 =  aux10*xn4
     b11 =  aux11*xn4
     b12 =  aux12*xn4
     b13 =  4.d0*(g2*aux1 + g1*aux4)
     b14 =  4.d0*(g2*aux2 + g1*aux5)
     b15 =  4.d0*(g2*aux3 + g1*aux6)
     b16 =  4.d0*(g3*aux4 + g2*aux7)
     b17 =  4.d0*(g3*aux5 + g2*aux8)
     b18 =  4.d0*(g3*aux6 + g2*aux9)
     b19 =  4.d0*(g1*aux7 + g3*aux1)
     b20 =  4.d0*(g1*aux8 + g3*aux2)
     b21 =  4.d0*(g1*aux9 + g3*aux3)
     b22 =  4.d0*(g4*aux1 + g1*aux10)
     b23 =  4.d0*(g4*aux2 + g1*aux11)
     b24 =  4.d0*(g4*aux3 + g1*aux12)
     b25 =  4.d0*(g4*aux4 + g2*aux10)
     b26 =  4.d0*(g4*aux5 + g2*aux11)
     b27 =  4.d0*(g4*aux6 + g2*aux12)
     b28 =  4.d0*(g4*aux7 + g3*aux10)
     b29 =  4.d0*(g4*aux8 + g3*aux11)
     b30 =  4.d0*(g4*aux9 + g3*aux12)

     e11 = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + &
          b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10) * vx6inv
     e22 = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + &
          b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10) * vx6inv
     e33 = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + &
          b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10) * vx6inv
     e12 = (b2*u1 + b1*v1 + b5*u2 + b4*v2  &
          + b8*u3 + b7*v3 + b11*u4 + b10*v4 &
          + b14*u5 + b13*v5 + b17*u6 + b16*v6 &
          + b20*u7 + b19*v7 + b23*u8 + b22*v8 &
          + b26*u9 + b25*v9 + b29*u10 + b28*v10) * vx6inv
     e23 = (b3*v1 + b2*w1 + b6*v2 + b5*w2  &
          + b9*v3 + b8*w3 + b12*v4 + b11*w4 &
          + b15*v5 + b14*w5 + b18*v6 + b17*w6 &
          + b21*v7 + b20*w7 + b24*v8 + b23*w8 &
          + b27*v9 + b26*w9 + b30*v10 + b29*w10) * vx6inv
     e13 = (b3*u1 + b1*w1 + b6*u2 + b4*w2  &
          + b9*u3 + b7*w3 + b12*u4 + b10*w4 &
          + b15*u5 + b13*w5 + b18*u6 + b16*w6 &
          + b21*u7 + b19*w7 + b24*u8 + b22*w8 &
          + b27*u9 + b25*w9 + b30*u10 + b28*w10) * vx6inv
     
     s11l(3,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
     s22l(3,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
     s33l(3,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
     s12l(3,i) = e12*ci(7,j)
     s23l(3,i) = e23*ci(8,j)
     s13l(3,i) = e13*ci(9,j)
     
     r1 =  r1 + s11l(3,i)*b1 + s12l(3,i)*b2 + s13l(3,i)*b3
     r2 =  r2 + s22l(3,i)*b2 + s12l(3,i)*b1 + s23l(3,i)*b3
     r3 =  r3 +s33l(3,i)*b3 + s23l(3,i)*b2 + s13l(3,i)*b1
     r4 =  r4 +s11l(3,i)*b4 + s12l(3,i)*b5 + s13l(3,i)*b6
     r5 =  r5 +s22l(3,i)*b5 + s12l(3,i)*b4 + s23l(3,i)*b6
     r6 =  r6 +s33l(3,i)*b6 + s23l(3,i)*b5 + s13l(3,i)*b4
     r7 =  r7 +s11l(3,i)*b7 + s12l(3,i)*b8 + s13l(3,i)*b9
     r8 =  r8 +s22l(3,i)*b8 + s12l(3,i)*b7 + s23l(3,i)*b9
     r9 =  r9 +s33l(3,i)*b9 + s23l(3,i)*b8 + s13l(3,i)*b7
     r10 = r10 +s11l(3,i)*b10 + s12l(3,i)*b11 + s13l(3,i)*b12
     r11 = r11 +s22l(3,i)*b11 + s12l(3,i)*b10 + s23l(3,i)*b12
     r12 = r12 +s33l(3,i)*b12 + s23l(3,i)*b11 + s13l(3,i)*b10
     r13 = r13 +s11l(3,i)*b13 + s12l(3,i)*b14 + s13l(3,i)*b15
     r14 = r14 +s22l(3,i)*b14 + s12l(3,i)*b13 + s23l(3,i)*b15
     r15 = r15 +s33l(3,i)*b15 + s23l(3,i)*b14 + s13l(3,i)*b13
     r16 = r16 +s11l(3,i)*b16 + s12l(3,i)*b17 + s13l(3,i)*b18
     r17 = r17 +s22l(3,i)*b17 + s12l(3,i)*b16 + s23l(3,i)*b18
     r18 = r18 +s33l(3,i)*b18 + s23l(3,i)*b17 + s13l(3,i)*b16
     r19 = r19 +s11l(3,i)*b19 + s12l(3,i)*b20 + s13l(3,i)*b21
     r20 = r20 +s22l(3,i)*b20 + s12l(3,i)*b19 + s23l(3,i)*b21
     r21 = r21 +s33l(3,i)*b21 + s23l(3,i)*b20 + s13l(3,i)*b19
     r22 = r22 +s11l(3,i)*b22 + s12l(3,i)*b23 + s13l(3,i)*b24
     r23 = r23 +s22l(3,i)*b23 + s12l(3,i)*b22 + s23l(3,i)*b24
     r24 = r24 +s33l(3,i)*b24 + s23l(3,i)*b23 + s13l(3,i)*b22
     r25 = r25 +s11l(3,i)*b25 + s12l(3,i)*b26 + s13l(3,i)*b27
     r26 = r26 +s22l(3,i)*b26 + s12l(3,i)*b25 + s23l(3,i)*b27
     r27 = r27 +s33l(3,i)*b27 + s23l(3,i)*b26 + s13l(3,i)*b25
     r28 = r28 +s11l(3,i)*b28 + s12l(3,i)*b29 + s13l(3,i)*b30
     r29 = r29 +s22l(3,i)*b29 + s12l(3,i)*b28 + s23l(3,i)*b30
     r30 = r30 +s33l(3,i)*b30 + s23l(3,i)*b29 + s13l(3,i)*b28
!--  Gauss Point 4

     g1 = 0.13819660d0
     g2 = 0.13819660d0
     g3 = 0.13819660d0
     g4 = 0.58541020d0

     xn1 = (4.d0*g1-1.d0)   ! derivative of shape function
     xn2 = (4.d0*g2-1.d0)   ! dN_i/dzeta_i
     xn3 = (4.d0*g3-1.d0)
     xn4 = (4.d0*g4-1.d0)
!        xN5 = 4.d0*g1*g2
!        xN6 = 4.d0*g2*g3
!        xN7 = 4.d0*g3*g1
!        xN8 = 4.d0*g1*g4
!        xN9 = 4.d0*g2*g4
!        xN10= 4.d0*g3*g4
     
     b1 = aux1*xn1
     b2 = aux2*xn1
     b3 = aux3*xn1
     b4 = aux4*xn2
     b5 = aux5*xn2
     b6 = aux6*xn2
     b7 = aux7*xn3
     b8 = aux8*xn3
     b9 = aux9*xn3
     b10 =  aux10*xn4
     b11 =  aux11*xn4
     b12 =  aux12*xn4
     b13 =  4.d0*(g2*aux1 + g1*aux4)
     b14 =  4.d0*(g2*aux2 + g1*aux5)
     b15 =  4.d0*(g2*aux3 + g1*aux6)
     b16 =  4.d0*(g3*aux4 + g2*aux7)
     b17 =  4.d0*(g3*aux5 + g2*aux8)
     b18 =  4.d0*(g3*aux6 + g2*aux9)
     b19 =  4.d0*(g1*aux7 + g3*aux1)
     b20 =  4.d0*(g1*aux8 + g3*aux2)
     b21 =  4.d0*(g1*aux9 + g3*aux3)
     b22 =  4.d0*(g4*aux1 + g1*aux10)
     b23 =  4.d0*(g4*aux2 + g1*aux11)
     b24 =  4.d0*(g4*aux3 + g1*aux12)
     b25 =  4.d0*(g4*aux4 + g2*aux10)
     b26 =  4.d0*(g4*aux5 + g2*aux11)
     b27 =  4.d0*(g4*aux6 + g2*aux12)
     b28 =  4.d0*(g4*aux7 + g3*aux10)
     b29 =  4.d0*(g4*aux8 + g3*aux11)
     b30 =  4.d0*(g4*aux9 + g3*aux12)
     
     e11 = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + &
          b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10) * vx6inv
     e22 = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + &
          b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10) * vx6inv
     e33 = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + &
          b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10) * vx6inv
     e12 = (b2*u1 + b1*v1 + b5*u2 + b4*v2  &
          + b8*u3 + b7*v3 + b11*u4 + b10*v4 &
          + b14*u5 + b13*v5 + b17*u6 + b16*v6 &
          + b20*u7 + b19*v7 + b23*u8 + b22*v8 &
          + b26*u9 + b25*v9 + b29*u10 + b28*v10) * vx6inv
     e23 = (b3*v1 + b2*w1 + b6*v2 + b5*w2  &
          + b9*v3 + b8*w3 + b12*v4 + b11*w4 &
          + b15*v5 + b14*w5 + b18*v6 + b17*w6 &
          + b21*v7 + b20*w7 + b24*v8 + b23*w8 &
          + b27*v9 + b26*w9 + b30*v10 + b29*w10) * vx6inv
     e13 = (b3*u1 + b1*w1 + b6*u2 + b4*w2  &
          + b9*u3 + b7*w3 + b12*u4 + b10*w4 &
          + b15*u5 + b13*w5 + b18*u6 + b16*w6 &
          + b21*u7 + b19*w7 + b24*u8 + b22*w8 &
          + b27*u9 + b25*w9 + b30*u10 + b28*w10) * vx6inv 
     
     s11l(4,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
     s22l(4,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
     s33l(4,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
     s12l(4,i) = e12*ci(7,j)
     s23l(4,i) = e23*ci(8,j)
     s13l(4,i) = e13*ci(9,j)
     
     r1 =  r1 + s11l(4,i)*b1 + s12l(4,i)*b2 + s13l(4,i)*b3
     r2 =  r2 + s22l(4,i)*b2 + s12l(4,i)*b1 + s23l(4,i)*b3
     r3 =  r3 + s33l(4,i)*b3 + s23l(4,i)*b2 + s13l(4,i)*b1
     r4 =  r4 + s11l(4,i)*b4 + s12l(4,i)*b5 + s13l(4,i)*b6
     r5 =  r5 + s22l(4,i)*b5 + s12l(4,i)*b4 + s23l(4,i)*b6
     r6 =  r6 + s33l(4,i)*b6 + s23l(4,i)*b5 + s13l(4,i)*b4
     r7 =  r7 + s11l(4,i)*b7 + s12l(4,i)*b8 + s13l(4,i)*b9
     r8 =  r8 + s22l(4,i)*b8 + s12l(4,i)*b7 + s23l(4,i)*b9
     r9 =  r9 + s33l(4,i)*b9 + s23l(4,i)*b8 + s13l(4,i)*b7
     r10 = r10 +s11l(4,i)*b10 + s12l(4,i)*b11 + s13l(4,i)*b12
     r11 = r11 +s22l(4,i)*b11 + s12l(4,i)*b10 + s23l(4,i)*b12
     r12 = r12 +s33l(4,i)*b12 + s23l(4,i)*b11 + s13l(4,i)*b10
     r13 = r13 +s11l(4,i)*b13 + s12l(4,i)*b14 + s13l(4,i)*b15
     r14 = r14 +s22l(4,i)*b14 + s12l(4,i)*b13 + s23l(4,i)*b15
     r15 = r15 +s33l(4,i)*b15 + s23l(4,i)*b14 + s13l(4,i)*b13
     r16 = r16 +s11l(4,i)*b16 + s12l(4,i)*b17 + s13l(4,i)*b18
     r17 = r17 +s22l(4,i)*b17 + s12l(4,i)*b16 + s23l(4,i)*b18
     r18 = r18 +s33l(4,i)*b18 + s23l(4,i)*b17 + s13l(4,i)*b16
     r19 = r19 +s11l(4,i)*b19 + s12l(4,i)*b20 + s13l(4,i)*b21
     r20 = r20 +s22l(4,i)*b20 + s12l(4,i)*b19 + s23l(4,i)*b21
     r21 = r21 +s33l(4,i)*b21 + s23l(4,i)*b20 + s13l(4,i)*b19
     r22 = r22 +s11l(4,i)*b22 + s12l(4,i)*b23 + s13l(4,i)*b24
     r23 = r23 +s22l(4,i)*b23 + s12l(4,i)*b22 + s23l(4,i)*b24
     r24 = r24 +s33l(4,i)*b24 + s23l(4,i)*b23 + s13l(4,i)*b22
     r25 = r25 +s11l(4,i)*b25 + s12l(4,i)*b26 + s13l(4,i)*b27
     r26 = r26 +s22l(4,i)*b26 + s12l(4,i)*b25 + s23l(4,i)*b27
     r27 = r27 +s33l(4,i)*b27 + s23l(4,i)*b26 + s13l(4,i)*b25
     r28 = r28 +s11l(4,i)*b28 + s12l(4,i)*b29 + s13l(4,i)*b30
     r29 = r29 +s22l(4,i)*b29 + s12l(4,i)*b28 + s23l(4,i)*b30
     r30 = r30 +s33l(4,i)*b30 + s23l(4,i)*b29 + s13l(4,i)*b28

! Wi (i.e. weight) for 4 guass integration is 1/4         

! Wi * 1/6 because the volume of a reference tetrahedra in 
! volume coordinates is 1/6

! ASSEMBLE THE INTERNAL FORCE VECTOR
!
! local node 1
     r_in(k1n1)  = r_in(k1n1)  - r1*0.04166666666666667d0
     r_in(k2n1)  = r_in(k2n1)  - r2*0.04166666666666667d0
     r_in(k3n1)  = r_in(k3n1)  - r3*0.04166666666666667d0
! local node 2
     r_in(k1n2)  = r_in(k1n2)  - r4*0.04166666666666667d0
     r_in(k2n2)  = r_in(k2n2)  - r5*0.04166666666666667d0
     r_in(k3n2)  = r_in(k3n2)  - r6*0.04166666666666667d0
! local node 3
     r_in(k1n3)  = r_in(k1n3)  - r7*0.04166666666666667d0
     r_in(k2n3)  = r_in(k2n3)  - r8*0.04166666666666667d0
     r_in(k3n3)  = r_in(k3n3)  - r9*0.04166666666666667d0
! local node 4
     r_in(k1n4)  = r_in(k1n4)  - r10*0.04166666666666667d0
     r_in(k2n4)  = r_in(k2n4)  - r11*0.04166666666666667d0
     r_in(k3n4)  = r_in(k3n4)  - r12*0.04166666666666667d0
! local node 5
     r_in(k1n5)  = r_in(k1n5)  - r13*0.04166666666666667d0
     r_in(k2n5)  = r_in(k2n5)  - r14*0.04166666666666667d0
     r_in(k3n5)  = r_in(k3n5)  - r15*0.04166666666666667d0
! local node 6
     r_in(k1n6)  = r_in(k1n6)  - r16*0.04166666666666667d0
     r_in(k2n6)  = r_in(k2n6)  - r17*0.04166666666666667d0
     r_in(k3n6)  = r_in(k3n6)  - r18*0.04166666666666667d0
! local node 7
     r_in(k1n7)  = r_in(k1n7)  - r19*0.04166666666666667d0
     r_in(k2n7)  = r_in(k2n7)  - r20*0.04166666666666667d0
     r_in(k3n7)  = r_in(k3n7)  - r21*0.04166666666666667d0
! local node 8
     r_in(k1n8)  = r_in(k1n8)  - r22*0.04166666666666667d0
     r_in(k2n8)  = r_in(k2n8)  - r23*0.04166666666666667d0
     r_in(k3n8)  = r_in(k3n8)  - r24*0.04166666666666667d0
! local node 9
     r_in(k1n9)  = r_in(k1n9)  - r25*0.04166666666666667d0
     r_in(k2n9)  = r_in(k2n9)  - r26*0.04166666666666667d0
     r_in(k3n9)  = r_in(k3n9)  - r27*0.04166666666666667d0
! local node 10
     r_in(k1n10) = r_in(k1n10) - r28*0.04166666666666667d0
     r_in(k2n10) = r_in(k2n10) - r29*0.04166666666666667d0
     r_in(k3n10) = r_in(k3n10) - r30*0.04166666666666667d0
  ENDDO
  
  RETURN
END SUBROUTINE v3d10