v3d10_B_Bar.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_b_bar(coor,matcstet,lmcstet,R_in,d,ci, &
     s11l,s22l,s33l,s12l,s23l,s13l,&
     numnp,nstart,nend,numcstet,numat_vol)
!-----Performs displacement based computations for 3D, 10-node 
!-----Quadratic Linear Elastic Tetrahedron with quadratic interpolation 
!-----functions. (linear strain tetrahedra)

!----- 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
!--   mat number for CSTet element
  INTEGER, DIMENSION(1:numcstet) :: matcstet
!--   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
!---- 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
  REAL*8 :: b4n1, b5n1, b6n1, b7n1, b8n1, b9n1
  REAL*8 :: b4n2, b5n2, b6n2, b7n2, b8n2, b9n2
  REAL*8 :: b4n3, b5n3, b6n3, b7n3, b8n3, b9n3
  REAL*8 :: b4n4, b5n4, b6n4, b7n4, b8n4, b9n4
  REAL*8 :: b4n5, b5n5, b6n5, b7n5, b8n5, b9n5
  REAL*8 :: b4n6, b5n6, b6n6, b7n6, b8n6, b9n6
  REAL*8 :: b4n7, b5n7, b6n7, b7n7, b8n7, b9n7
  REAL*8 :: b4n8, b5n8, b6n8, b7n8, b8n8, b9n8
  REAL*8 :: b4n9, b5n9, b6n9, b7n9, b8n9, b9n9
  REAL*8 :: b4n10, b5n10, b6n10, b7n10, b8n10, b9n10
!--   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,k
  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
!--  Coordinate subtractions: These are to speed up B calculation
  REAL*8 :: x12, x13, y12, y13, z12, z13
!--  
  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

  REAL*8 :: onethird = 1./3.

  REAL*8,DIMENSION(1:4,1:4) :: gaussintpt = reshape( &
       (/0.58541020d0,0.13819660d0,0.13819660d0,0.13819660d0, &
         0.13819660d0,0.58541020d0,0.13819660d0,0.13819660d0, &
         0.13819660d0,0.13819660d0,0.58541020d0,0.13819660d0, &
         0.13819660d0,0.13819660d0,0.13819660d0,0.58541020d0/),(/4,4/) )

  integer :: igpt

  REAL*8, DIMENSION(1:10) :: bbarx, bbary, bbarz
  REAL*8, DIMENSION(1:3,1:3) :: jac
  REAL*8, DIMENSION(1:3) :: dx,dy,dz
  REAL*8 :: jacinv

! Derivatives of the shape functions evaluated at the center of the tet, 1/4,1/4,1/4

  REAL*8, DIMENSION(1:10) :: n_ksi = (/0.,0.,0.,0.,0.,1.,-1.,-1.,1.,0./)
  REAL*8, DIMENSION(1:10) :: n_eta = (/0.,0.,0.,0.,-1.,1.,0.,-1.,0.,1./)
  REAL*8, DIMENSION(1:10) :: n_zet = (/0.,0.,0.,0.,-1.,0.,-1.,0.,1.,1./)


  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) 
     x5 = coor(1,n5)
     x6 = coor(1,n6)
     x7 = coor(1,n7)
     x8 = coor(1,n8)
     x9  = coor(1,n9)
     x10 = coor(1,n10)

     y1 = coor(2,n1)
     y2 = coor(2,n2)
     y3 = coor(2,n3)
     y4 = coor(2,n4)
     y5 = coor(2,n5)
     y6 = coor(2,n6)
     y7 = coor(2,n7)
     y8 = coor(2,n8)
     y9 = coor(2,n9)
     y10 = coor(2,n10)

     z1 = coor(3,n1)
     z2 = coor(3,n2)
     z3 = coor(3,n3)
     z4 = coor(3,n4)
     z5 = coor(3,n5)
     z6 = coor(3,n6)
     z7 = coor(3,n7)
     z8 = coor(3,n8)
     z9 = coor(3,n9)
     z10 = coor(3,n10)

     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

     x12 = x1 - x2          ! not used in vol. calc
     x13 = x1 - x3          ! not used in vol. calc
     y12 = y1 - y2          ! not used in vol. calc
     y13 = y1 - y3          ! not used in vol. calc
     z12 = z1 - z2          ! not used in vol. calc
     z13 = z1 - z3          ! not used in vol. calc

     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)


!!$     dx(1) =  x6 - x7 - x8 + x9  ! dx(1) + N_ksi(i)*x(i)
!!$     dx(2) = -x5 + x6 - x8 + x10  ! N_eta(i)*x(i)
!!$     dx(3) = -x5 - x7 + x9 + x10  ! N_zet(i)*x(i)
!!$     dy(1) =  y6 - y7 - y8 + y9 ! N_ksi(i)*y(i)
!!$     dy(2) = -y5 + y6 - y8 + y10 ! N_eta(i)*y(i)
!!$     dy(3) = -y5 - y7 + y9 + y10 ! N_zet(i)*y(i)
!!$     dz(1) =  z6 - z7 - z8 + z9 ! N_ksi(i)*z(i)
!!$     dz(2) = -z5 + z6 - z8 + z10 ! N_eta(i)*z(i)
!!$     dz(3) = -z5 - z7 + z9 + z10 ! N_zet(i)*z(i)
!!$
!!$     Jac(1,1) = dy(2)*dz(3) - dz(2)*dy(3)
!!$     Jac(1,2) = dz(1)*dy(3) - dy(1)*dz(3)
!!$     Jac(1,3) = dy(1)*dz(2) - dz(1)*dy(2)
!!$     Jac(2,1) = dz(2)*dx(3) - dx(2)*dz(3)
!!$     Jac(2,2) = dx(1)*dz(3) - dz(1)*dx(3)
!!$     Jac(2,3) = dz(1)*dx(2) - dx(1)*dz(2)
!!$     Jac(3,1) = dx(2)*dy(3) - dy(2)*dx(3)
!!$     Jac(3,2) = dy(1)*dx(3) - dx(1)*dy(3)
!!$     Jac(3,3) = dx(1)*dy(2) - dy(1)*dx(2)

!!$     JacInv =  1./( (dx(1)*dy(2)*dz(3)) + (dy(1)*dz(2)*dx(3)) + (dz(1)*dx(2)*dy(3)) & 
!!$          - (dz(1)*dy(2)*dx(3)) - (dx(1)*dz(2)*dy(3)) - (dy(1)*dx(2)*dz(3)))
!!$
!!$     PRINT*,JacInv,Vx6inv

!!$     JacInv = 1.0
!!$
!!$
!!$     DO k = 1, 10
!!$        BBarX(k) = JacInv*(Jac(1,1)*N_ksi(k) + Jac(1,2)*N_eta(k) + Jac(1,3)*N_zet(k))
!!$        BBarY(k) = JacInv*(Jac(2,1)*N_ksi(k) + Jac(2,2)*N_eta(k) + Jac(2,3)*N_zet(k))
!!$        BBarZ(k) = JacInv*(Jac(3,1)*N_ksi(k) + Jac(3,2)*N_eta(k) + Jac(3,3)*N_zet(k))
!!$     ENDDO

     bbarx(1) = 0. ! aux1*N_ksi(1)
     bbary(1) = 0. !aux2*N_eta(1)
     bbarz(1) = 0. !aux3*N_zet(1)
     bbarx(2) = 0. !aux4*N_ksi(2)
     bbary(2) = 0. !aux5*N_eta(2)
     bbarz(2) = 0. !aux6*N_zet(2)
     bbarx(3) = 0. !aux4*N_ksi(3)
     bbary(3) = 0. !aux5*N_eta(3)
     bbarz(3) = 0. !aux6*N_zet(3)
     bbarx(4) = 0. !aux7*N_ksi(4)
     bbary(4) = 0. !aux8*N_eta(4)
     bbarz(4) = 0. !aux9*N_zet(4)

     bbarx(5) =  aux1 + aux4
     bbary(5) =  aux2 + aux5
     bbarz(5) =  aux3 + aux6

     bbarx(6) =  aux4 + aux7
     bbary(6) =  aux5 + aux8
     bbarz(6) =  aux6 + aux9
        
     bbarx(7) =  aux7 + aux1
     bbary(7) =  aux8 + aux2
     bbarz(7) =  aux9 + aux3
     
     bbarx(8) =  aux1 + aux10
     bbary(8) =  aux2 + aux11
     bbarz(8) =  aux3 + aux12
     
     bbarx(9) =  aux4 + aux10
     bbary(9) =  aux5 + aux11
     bbarz(9) =  aux6 + aux12
     
     bbarx(10) =  aux7 + aux10
     bbary(10) =  aux8 + aux11
     bbarz(10) =  aux9 + aux12

     r1 = 0.
     r2 = 0.
     r3 = 0.
     r4 = 0.
     r5 = 0.
     r6 = 0.
     r7 = 0.
     r8 = 0.
     r9 = 0.
     r10 = 0.
     r11 = 0.
     r12 = 0.
     r13 = 0.
     r14 = 0.
     r15 = 0.
     r16 = 0.
     r17 = 0.
     r18 = 0.
     r19 = 0.
     r20 = 0.
     r21 = 0.
     r22 = 0.
     r23 = 0.
     r24 = 0.
     r25 = 0.
     r26 = 0.
     r27 = 0.
     r28 = 0.
     r29 = 0.
     r30 = 0.


     DO igpt = 1, 4

        g1 = gaussintpt(igpt,1)
        g2 = gaussintpt(igpt,2)
        g3 = gaussintpt(igpt,3)
        g4 = gaussintpt(igpt,4)


        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)

!!$        BBarX(1) = B1
!!$        BBarY(1) = B2
!!$        BBarZ(1) = B3
!!$        BBarX(2) = B4
!!$        BBarY(2) = B5
!!$        BBarZ(2) = B6
!!$        BBarX(3) = B7
!!$        BBarY(3) = B8
!!$        BBarZ(3) = B9
!!$        BBarX(4) = B10
!!$        BBarY(4) = B11
!!$        BBarZ(4) = B12
!!$        BBarX(5) = B13
!!$        BBarY(5) = B14
!!$        BBarZ(5) = B15
!!$        BBarX(6) = B16
!!$        BBarY(6) = B17
!!$        BBarZ(6) = B18
!!$        BBarX(7) = B19
!!$        BBarY(7) = B20
!!$        BBarZ(7) = B21
!!$        BBarX(8) = B22
!!$        BBarY(8) = B23
!!$        BBarZ(8) = B24
!!$        BBarX(9) = B25
!!$        BBarY(9) = B26
!!$        BBarZ(9) = B27
!!$        BBarX(10) = B28
!!$        BBarY(10) = B29
!!$        BBarZ(10) = B30
     
        b4n1 = ( bbarx(1) - b1 )*onethird
        b5n1 = b1 + b4n1
        b6n1 = ( bbary(1) - b2 )*onethird
        b7n1 = b2 + b6n1
        b8n1 = ( bbarz(1) - b3 )*onethird
        b9n1 = b3 + b8n1
        
        b4n2 = ( bbarx(2) - b4 )*onethird
        b5n2 = b4 + b4n2
        b6n2 = ( bbary(2) - b5 )*onethird
        b7n2 = b5 + b6n2
        b8n2 = ( bbarz(2) - b6 )*onethird
        b9n2 = b6 + b8n2
        
        b4n3 = ( bbarx(3) - b7 )*onethird
        b5n3 = b7 + b4n3
        b6n3 = ( bbary(3) - b8 )*onethird
        b7n3 = b8 + b6n3
        b8n3 = ( bbarz(3) - b9 )*onethird
        b9n3 = b9 + b8n3
        
        b4n4 = ( bbarx(4) - b10 )*onethird
        b5n4 = b10 + b4n4
        b6n4 = ( bbary(4) - b11 )*onethird
        b7n4 = b11 + b6n4
        b8n4 = ( bbarz(4) - b12 )*onethird
        b9n4 = b12 + b8n4
        
        b4n5 = ( bbarx(5) - b13 )*onethird
        b5n5 = b13 + b4n5
        b6n5 = ( bbary(5) - b14 )*onethird
        b7n5 = b14 + b6n5
        b8n5 = ( bbarz(5) - b15 )*onethird
        b9n5 = b15 + b8n5   
        
        b4n6 = ( bbarx(6) - b16 )*onethird
        b5n6 = b16 + b4n6
        b6n6 = ( bbary(6) - b17 )*onethird
        b7n6 = b17 + b6n6
        b8n6 = ( bbarz(6) - b18 )*onethird
        b9n6 = b18 + b8n6
        
        b4n7 = ( bbarx(7) - b19 )*onethird
        b5n7 = b19 + b4n7
        b6n7 = ( bbary(7) - b20 )*onethird
        b7n7 = b20 + b6n7
        b8n7 = ( bbarz(7) - b21 )*onethird
        b9n7 = b21 + b8n7
        
        b4n8 = ( bbarx(8) - b22 )*onethird
        b5n8 = b22 + b4n8
        b6n8 = ( bbary(8) - b23 )*onethird
        b7n8 = b23 + b6n8
        b8n8 = ( bbarz(8) - b24 )*onethird
        b9n8 = b24 + b8n8
        
        b4n9 = ( bbarx(9) - b25 )*onethird
        b5n9 = b25 + b4n9
        b6n9 = ( bbary(9) - b26 )*onethird
        b7n9 = b26 + b6n9
        b8n9 = ( bbarz(9) - b27 )*onethird
        b9n9 = b27 + b8n9
        
        b4n10 = ( bbarx(10) - b28 )*onethird
        b5n10 = b28 + b4n10
        b6n10 = ( bbary(10) - b29 )*onethird
        b7n10 = b29 + b6n10
        b8n10 = ( bbarz(10) - b30 )*onethird
        b9n10 = b30 + b8n10
        
        e11 = (b5n1*u1 + b6n1*v1 + b8n1*w1 + b5n2*u2 + b6n2*v2 + b8n2*w2 + &
             b5n3*u3 + b6n3*v3 + b8n3*w3 + b5n4*u4 + b6n4*v4 + b8n4*w4 + &
             b5n5*u5 + b6n5*v5 + b8n5*w5 + b5n6*u6 + b6n6*v6 + b8n6*w6 + &
             b5n7*u7 + b6n7*v7 + b8n7*w7 + b5n8*u8 + b6n8*v8 + b8n8*w8 + &
             b5n9*u9 + b6n9*v9 + b8n9*w9 + b5n10*u10 + b6n10*v10 + b8n10*w10) * vx6inv
        
        e22 = (b4n1*u1 + b7n1*v1 + b8n1*w1 + b4n2*u2 + b7n2*v2 + b8n2*w2 +  &
             b4n3*u3 + b7n3*v3 + b8n3*w3 + b4n4*u4 + b7n4*v4 + b8n4*w4 +  &
             b4n5*u5 + b7n5*v5 + b8n5*w5 + b4n6*u6 + b7n6*v6 + b8n6*w6 +  &
             b4n7*u7 + b7n7*v7 + b8n7*w7 + b4n8*u8 + b7n8*v8 + b8n8*w8 +  &
             b4n9*u9 + b7n9*v9 + b8n9*w9 + b4n10*u10 + b7n10*v10 + b8n10*w10) * vx6inv
        
        e33 = (b4n1*u1 + b6n1*v1 + b9n1*w1 + b4n2*u2 + b6n2*v2 + b9n2*w2 +  &
             b4n3*u3 + b6n3*v3 + b9n3*w3 + b4n4*u4 + b6n4*v4 + b9n4*w4 +  &
             b4n5*u5 + b6n5*v5 + b9n5*w5 + b4n6*u6 + b6n6*v6 + b9n6*w6 +  &
             b4n7*u7 + b6n7*v7 + b9n7*w7 + b4n8*u8 + b6n8*v8 + b9n8*w8 +  &
             b4n9*u9 + b6n9*v9 + b9n9*w9 + b4n10*u10 + b6n10*v10 + b9n10*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(igpt,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
        s22l(igpt,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
        s33l(igpt,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
        s12l(igpt,i) = e12*ci(7,j)
        s23l(igpt,i) = e23*ci(8,j)
        s13l(igpt,i) = e13*ci(9,j)
     
        r1 = r1 + s11l(igpt,i)*b5n1 + s22l(igpt,i)*b4n1 + s33l(igpt,i)*b4n1 + s12l(igpt,i)*b2 + s13l(igpt,i)*b3
        r2 = r2 +s11l(igpt,i)*b6n1 + s22l(igpt,i)*b7n1 + s33l(igpt,i)*b6n1 + s12l(igpt,i)*b1 + s23l(igpt,i)*b3
        r3 = r3 +s11l(igpt,i)*b8n1 + s22l(igpt,i)*b8n1 + s33l(igpt,i)*b9n1 + s23l(igpt,i)*b2 + s13l(igpt,i)*b1
        
        r4 = r4 +s11l(igpt,i)*b5n2 + s22l(igpt,i)*b4n2 + s33l(igpt,i)*b4n2 + s12l(igpt,i)*b5 + s13l(igpt,i)*b6
        r5 = r5 +s11l(igpt,i)*b6n2 + s22l(igpt,i)*b7n2 + s33l(igpt,i)*b6n2 + s12l(igpt,i)*b4 + s23l(igpt,i)*b6
        r6 = r6 +s11l(igpt,i)*b8n2 + s22l(igpt,i)*b8n2 + s33l(igpt,i)*b9n2 + s23l(igpt,i)*b5 + s13l(igpt,i)*b4
        
        r7 = r7 +s11l(igpt,i)*b5n3 + s22l(igpt,i)*b4n3 + s33l(igpt,i)*b4n3 + s12l(igpt,i)*b8 + s13l(igpt,i)*b9
        r8 = r8 +s11l(igpt,i)*b6n3 + s22l(igpt,i)*b7n3 + s33l(igpt,i)*b6n3 + s12l(igpt,i)*b7 + s23l(igpt,i)*b9
        r9 = r9 +s11l(igpt,i)*b8n3 + s22l(igpt,i)*b8n3 + s33l(igpt,i)*b9n3 + s23l(igpt,i)*b8 + s13l(igpt,i)*b7
        
        r10 = r10 +s11l(igpt,i)*b5n4 + s22l(igpt,i)*b4n4 + s33l(igpt,i)*b4n4 + s12l(igpt,i)*b11 + s13l(igpt,i)*b12
        r11 = r11 +s11l(igpt,i)*b6n4 + s22l(igpt,i)*b7n4 + s33l(igpt,i)*b6n4 + s12l(igpt,i)*b10 + s23l(igpt,i)*b12
        r12 = r12 +s11l(igpt,i)*b8n4 + s22l(igpt,i)*b8n4 + s33l(igpt,i)*b9n4 + s23l(igpt,i)*b11 + s13l(igpt,i)*b10
        
        r13 = r13 +s11l(igpt,i)*b5n5 + s22l(igpt,i)*b4n5 + s33l(igpt,i)*b4n5 + s12l(igpt,i)*b14 + s13l(igpt,i)*b15
        r14 = r14 +s11l(igpt,i)*b6n5 + s22l(igpt,i)*b7n5 + s33l(igpt,i)*b6n5 + s12l(igpt,i)*b13 + s23l(igpt,i)*b15
        r15 = r15 +s11l(igpt,i)*b8n5 + s22l(igpt,i)*b8n5 + s33l(igpt,i)*b9n5 + s23l(igpt,i)*b14 + s13l(igpt,i)*b13
        
        r16 = r16 +s11l(igpt,i)*b5n6 + s22l(igpt,i)*b4n6 + s33l(igpt,i)*b4n6 + s12l(igpt,i)*b17 + s13l(igpt,i)*b18
        r17 = r17 +s11l(igpt,i)*b6n6 + s22l(igpt,i)*b7n6 + s33l(igpt,i)*b6n6 + s12l(igpt,i)*b16 + s23l(igpt,i)*b18
        r18 = r18 +s11l(igpt,i)*b8n6 + s22l(igpt,i)*b8n6 + s33l(igpt,i)*b9n6 + s23l(igpt,i)*b17 + s13l(igpt,i)*b16
        
        r19 = r19 +s11l(igpt,i)*b5n7 + s22l(igpt,i)*b4n7 + s33l(igpt,i)*b4n7 + s12l(igpt,i)*b20 + s13l(igpt,i)*b21
        r20 = r20 +s11l(igpt,i)*b6n7 + s22l(igpt,i)*b7n7 + s33l(igpt,i)*b6n7 + s12l(igpt,i)*b19 + s23l(igpt,i)*b21
        r21 = r21 +s11l(igpt,i)*b8n7 + s22l(igpt,i)*b8n7 + s33l(igpt,i)*b9n7 + s23l(igpt,i)*b20 + s13l(igpt,i)*b19
        
        r22 = r22 +s11l(igpt,i)*b5n8 + s22l(igpt,i)*b4n8 + s33l(igpt,i)*b4n8 + s12l(igpt,i)*b23 + s13l(igpt,i)*b24
        r23 = r23 +s11l(igpt,i)*b6n8 + s22l(igpt,i)*b7n8 + s33l(igpt,i)*b6n8 + s12l(igpt,i)*b22 + s23l(igpt,i)*b24
        r24 = r24 +s11l(igpt,i)*b8n8 + s22l(igpt,i)*b8n8 + s33l(igpt,i)*b9n8 + s23l(igpt,i)*b23 + s13l(igpt,i)*b22
        
        r25 = r25 +s11l(igpt,i)*b5n9 + s22l(igpt,i)*b4n9 + s33l(igpt,i)*b4n9 + s12l(igpt,i)*b26 + s13l(igpt,i)*b27
        r26 = r26 +s11l(igpt,i)*b6n9 + s22l(igpt,i)*b7n9 + s33l(igpt,i)*b6n9 + s12l(igpt,i)*b25 + s23l(igpt,i)*b27
        r27 = r27 +s11l(igpt,i)*b8n9 + s22l(igpt,i)*b8n9 + s33l(igpt,i)*b9n9 + s23l(igpt,i)*b26 + s13l(igpt,i)*b25
        
        r28 = r28 +s11l(igpt,i)*b5n10 + s22l(igpt,i)*b4n10 + s33l(igpt,i)*b4n10 + s12l(igpt,i)*b29 + s13l(igpt,i)*b30
        r29 = r29 +s11l(igpt,i)*b6n10 + s22l(igpt,i)*b7n10 + s33l(igpt,i)*b6n10 + s12l(igpt,i)*b28 + s23l(igpt,i)*b30
        r30 = r30 +s11l(igpt,i)*b8n10 + s22l(igpt,i)*b8n10 + s33l(igpt,i)*b9n10 + s23l(igpt,i)*b29 + s13l(igpt,i)*b28
        
     ENDDO

! 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_b_bar