v3d4n_NeoHookeanCompress.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 v3d4n_neohookeancompress(numnp, numel, coor, disp, nodes, Rnet,&
     numelneigh, elconn, alpha, ahat,nummatvol, &
     xmu,xlambda)
  
!     nprocs,TotNumNdComm,TotNumNeighProcs,NeighProcList,NumNdComm,neigh_lst,&
!     MPI_STATUS_SIZE,MPI_COMM_ROCFRAC,MPI_DOUBLE_PRECISION, &
!     ReqRcv,ReqSnd,StatRcv, StatSnd)


!________________________________________________________________________
!
!  V3D4 - Performs displacement based computations for Volumetric 3D 
!         4-node tetrahedra large deformation (i.e. nonlinear) elastic 
!         elements with linear interpolation functions. (constant strain 
!         tetrahedra). Returns the internal force vector R_in.
!
!  DATE: 04.2000                  AUTHOR: SCOT BREITENFELD
!
! Source:
!    "Nonlinear continuum mechanics for finite element analysis"
!     Javier Bonet and Richard D. Wood
!     ISBN 0-521-57272-X

!________________________________________________________________________

  USE rocstar_rocfraccomm 

  IMPLICIT NONE
  
  integer :: nummatvol
  INTEGER :: numnp, numel
  REAL*8,dimension(1:NumMatVol) :: xmu, xlambda
  REAL*8, DIMENSION(1:numnp) :: ahat
  REAL*8, DIMENSION(1:3*numnp) :: disp, rnet
  REAL*8, DIMENSION(1:3,1:numnp) :: coor
  INTEGER, DIMENSION(1:numnp) ::numelneigh
  INTEGER, DIMENSION(1:numnp,1:40) :: elconn  ! fix 40 should not be hard coded
  INTEGER,DIMENSION(1:4,1:numel) :: nodes
  REAL*8, DIMENSION(1:4,1:numel) :: alpha
  REAL*8, DIMENSION(1:numnp) :: s11n, s22n, s33n, s12n, s13n,s23n

  INTEGER :: i,k
!--  x, y and z displacements of nodes
  REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
!--   coordinate holding variable
  REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
!--  Coordinate subtractions

  REAL*8 ::  z12, z13
  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
  REAL*8 :: x12, x13, y12, y13

!--   6*volume, inverse(6*volume),  and the volume      
  REAL*8 :: vx6, vx6inv, vol
!--   spacial derivatives
  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
!--   partial derivatives of the displacement 
  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
!--   strains
  REAL*8 :: volnd,vainv
! -- node numbers
  INTEGER :: n1,n2,n3,n4
! -- dummy and counters
  INTEGER :: j,nstart,nend
  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
  INTEGER :: k3n1,k3n2,k3n3,k3n4
! -- 
  REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
  REAL*8 :: j1,j2,j2inv
  REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33
  INTEGER :: jj,ielnum
  REAL*8 :: volnd0inv,volel0
  REAL*8, DIMENSION(1:1) :: s11, s22, s33, s12, s23, s13
  INTEGER :: k3i,k2i,k1i
  INTEGER :: ix

! /* For each Node */

  DO i = 1, numnp ! for each node

! /* Undeformed Volume of node */

     volnd0inv = 1.d0/ahat(i)

! /* Initialize to 0, Va, Va0, Fa */

     f11 = 0.d0
     f22 = 0.d0
     f33 = 0.d0
     f12 = 0.d0
     f13 = 0.d0
     f21 = 0.d0
     f23 = 0.d0
     f31 = 0.d0
     f32 = 0.d0
     volnd = 0.d0
       
     DO j = 1, numelneigh(i)
        
        ielnum = elconn(i,j)

        n1 = nodes(1,ielnum)
        n2 = nodes(2,ielnum)
        n3 = nodes(3,ielnum)
        n4 = nodes(4,ielnum)

        IF(i.EQ.n1) THEN
           ix = 1
        ELSE IF(i.EQ.n2)THEN
           ix = 2
        ELSE IF(i.EQ.n3)THEN
           ix = 3
        ELSE IF(i.EQ.n4)THEN
           ix = 4
        ENDIF

        k3n1 = 3*n1
        k3n2 = 3*n2
        k3n3 = 3*n3
        k3n4 = 3*n4
            
        k2n1 = k3n1 - 1
        k2n2 = k3n2 - 1
        k2n3 = k3n3 - 1
        k2n4 = k3n4 - 1

        k1n1 = k3n1 - 2
        k1n2 = k3n2 - 2
        k1n3 = k3n3 - 2
        k1n4 = k3n4 - 2     
                         ! k#n# dummy variables replaces:
        u1 = disp(k1n1)           ! (3*n1-2)
        u2 = disp(k1n2)           ! (3*n2-2)
        u3 = disp(k1n3)           ! (3*n3-2)
        u4 = disp(k1n4)           ! (3*n4-2)
        v1 = disp(k2n1)           ! (3*n1-1)
        v2 = disp(k2n2)           ! (3*n2-1)
        v3 = disp(k2n3)           ! (3*n3-1)
        v4 = disp(k2n4)           ! (3*n4-1)
        w1 = disp(k3n1)           ! (3*n1)
        w2 = disp(k3n2)           ! (3*n2)
        w3 = disp(k3n3)           ! (3*n3)
        w4 = disp(k3n4)           ! (3*n4)
            
        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)

        x12 = x1 - x2 ! not used in vol. calc
        x13 = x1 - x3 ! not used in vol. calc
        x14 = x1 - x4
        x24 = x2 - x4
        x34 = x3 - x4
        y12 = y1 - y2 ! not used in vol. calc
        y13 = y1 - y3 ! not used in vol. calc
        y14 = y1 - y4
        y24 = y2 - y4
        y34 = y3 - y4
        z12 = z1 - z2 ! not used in vol. calc
        z13 = z1 - z3 ! not used in vol. calc
        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 )

        vx6inv = 1.d0/vx6

        volel0 = vx6/6.d0 ! undeformed volume of element (Ve_O)
            

! See the maple worksheet 'V3D4.mws' for the derivation of Vx6


! See the maple worksheet 'V3D4.mws' for the derivation of [B]
! Compute the Shape functions
! NOTE: Factored for a more equivalent/compact form then maple's

        b1  = (y34*z24 - y24*z34) * vx6inv
        b2  = (z34*x24 - z24*x34) * vx6inv
        b3  = (x34*y24 - x24*y34) * vx6inv
        b4  = (y13*z14 - y14*z13) * vx6inv
        b5  = (z13*x14 - z14*x13) * vx6inv
        b6  = (x13*y14 - x14*y13) * vx6inv
        b7  = (y14*z12 - y12*z14) * vx6inv
        b8  = (z14*x12 - z12*x14) * vx6inv
        b9  = (x14*y12 - x12*y14) * vx6inv
        b10 = (y12*z13 - y13*z12) * vx6inv
        b11 = (z12*x13 - z13*x12) * vx6inv
        b12 = (x12*y13 - x13*y12) * vx6inv
!
! Calculate displacement gradient (H)
!
        dudx = b1*u1 + b4*u2 + b7*u3 + b10*u4
        dvdy = b2*v1 + b5*v2 + b8*v3 + b11*v4
        dwdz = b3*w1 + b6*w2 + b9*w3 + b12*w4
        dudy = b2*u1 + b5*u2 + b8*u3 + b11*u4
        dvdx = b1*v1 + b4*v2 + b7*v3 + b10*v4
        dvdz = b3*v1 + b6*v2 + b9*v3 + b12*v4
        dwdy = b2*w1 + b5*w2 + b8*w3 + b11*w4
        dudz = b3*u1 + b6*u2 + b9*u3 + b12*u4
        dwdx = b1*w1 + b4*w2 + b7*w3 + b10*w4
! 
! deformation gradients F
!
        f11 = f11 + alpha(ix,ielnum)*volel0*(1.d0 + dudx)
        f12 = f12 + alpha(ix,ielnum)*volel0*dudy
        f13 = f13 + alpha(ix,ielnum)*volel0*dudz
        f21 = f21 + alpha(ix,ielnum)*volel0*dvdx
        f22 = f22 + alpha(ix,ielnum)*volel0*(1.d0 + dvdy)
        f23 = f23 + alpha(ix,ielnum)*volel0*dvdz
        f31 = f31 + alpha(ix,ielnum)*volel0*dwdx
        f32 = f32 + alpha(ix,ielnum)*volel0*dwdy
        f33 = f33 + alpha(ix,ielnum)*volel0*(1.d0 + dwdz)

     ENDDO

! 2nd "for each node" in box starts here


  ! first bullet in box
     
     f11 = volnd0inv*f11
     f22 = volnd0inv*f22
     f33 = volnd0inv*f33
     f12 = volnd0inv*f12
     f21 = volnd0inv*f21
     f23 = volnd0inv*f23
     f32 = volnd0inv*f32
     f13 = volnd0inv*f13
     f31 = volnd0inv*f31

     j1 = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)
     IF(j1.LE.0.d0) THEN
        WRITE(*,*)'area has become zero for element',i
        stop
     ENDIF
     j2 = j1*j1

!          T
!     C = F  F = right Cauchy-Green deformation tensor
!
!     Divided by third invariant squared

     j2inv = 1.d0/j2
     
     c11 = (f11*f11+f21*f21+f31*f31)*j2inv
     c12 = (f11*f12+f21*f22+f31*f32)*j2inv
     c13 = (f11*f13+f21*f23+f31*f33)*j2inv
     c21 = c12
     c22 = (f12*f12+f22*f22+f32*f32)*j2inv
     c23 = (f12*f13+f22*f23+f32*f33)*j2inv
     c31 = c13
     c32 = c23
     c33 = (f13*f13+f23*f23+f33*f33)*j2inv       
     
!
! Second Piola-Kirchoff tensor
! Eq. (5.28), pg. 124

     s11n(i) = xmu(1)*(1.d0-(c22*c33-c23*c32)) + xlambda(1)*log(j1)*(c22*c33-c23*c32)
     s22n(i) = xmu(1)*(1.d0-(c11*c33-c31*c13)) + xlambda(1)*log(j1)*(c11*c33-c31*c13)
     s33n(i) = xmu(1)*(1.d0-(c11*c22-c12*c21)) + xlambda(1)*log(j1)*(c11*c22-c12*c21)
     s12n(i) = (-c12*c33+c13*c32)*(xmu(1) - xlambda(1)*log(j1))
     s23n(i) = (-c11*c23+c13*c21)*(xmu(1) - xlambda(1)*log(j1))
     s13n(i) = (c12*c23-c13*c22)*(xmu(1) - xlambda(1)*log(j1))
     
  ENDDO

! For each node

  DO i = 1, numnp

     k3i = 3*i
     k2i = 3*i-1
     k1i = 3*i-2

! For each element

     DO j = 1, numelneigh(i)
        
        ielnum = elconn(i,j)
        
        n1 = nodes(1,ielnum)
        n2 = nodes(2,ielnum)
        n3 = nodes(3,ielnum)
        n4 = nodes(4,ielnum)  
        k3n1 = 3*n1
        k3n2 = 3*n2
        k3n3 = 3*n3
        k3n4 = 3*n4
        
        k2n1 = k3n1 - 1
        k2n2 = k3n2 - 1
        k2n3 = k3n3 - 1
        k2n4 = k3n4 - 1
        
        k1n1 = k3n1 - 2
        k1n2 = k3n2 - 2
        k1n3 = k3n3 - 2
        k1n4 = k3n4 - 2     
        ! k#n# dummy variables replaces:
        u1 = disp(k1n1)           ! (3*n1-2)
        u2 = disp(k1n2)           ! (3*n2-2)
        u3 = disp(k1n3)           ! (3*n3-2)
        u4 = disp(k1n4)           ! (3*n4-2)
        v1 = disp(k2n1)           ! (3*n1-1)
        v2 = disp(k2n2)           ! (3*n2-1)
        v3 = disp(k2n3)           ! (3*n3-1)
        v4 = disp(k2n4)           ! (3*n4-1)
        w1 = disp(k3n1)           ! (3*n1)
        w2 = disp(k3n2)           ! (3*n2)
        w3 = disp(k3n3)           ! (3*n3)
        w4 = disp(k3n4)           ! (3*n4)
        
        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)
        
        x12 = x1 - x2 ! not used in vol. calc
        x13 = x1 - x3 ! not used in vol. calc
        x14 = x1 - x4
        x24 = x2 - x4
        x34 = x3 - x4
        y12 = y1 - y2 ! not used in vol. calc
        y13 = y1 - y3 ! not used in vol. calc
        y14 = y1 - y4
        y24 = y2 - y4
        y34 = y3 - y4
        z12 = z1 - z2 ! not used in vol. calc
        z13 = z1 - z3 ! not used in vol. calc
        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 )
        volel0 = vx6/6.d0

! See the maple worksheet 'V3D4.mws' for the derivation of Vx6

        vx6inv = 1.d0 / vx6

! See the maple worksheet 'V3D4.mws' for the derivation of [B]! Compute the Shape functions
! NOTE: Factored for a more equivalent/compact form then maple's

        b1  = (y34*z24 - y24*z34) * vx6inv
        b2  = (z34*x24 - z24*x34) * vx6inv
        b3  = (x34*y24 - x24*y34) * vx6inv
        b4  = (y13*z14 - y14*z13) * vx6inv
        b5  = (z13*x14 - z14*x13) * vx6inv
        b6  = (x13*y14 - x14*y13) * vx6inv
        b7  = (y14*z12 - y12*z14) * vx6inv
        b8  = (z14*x12 - z12*x14) * vx6inv
        b9  = (x14*y12 - x12*y14) * vx6inv
        b10 = (y12*z13 - y13*z12) * vx6inv
        b11 = (z12*x13 - z13*x12) * vx6inv
        b12 = (x12*y13 - x13*y12) * vx6inv
!
! Calculate displacement gradient (H)
!
        dudx = b1*u1 + b4*u2 + b7*u3 + b10*u4
        dvdy = b2*v1 + b5*v2 + b8*v3 + b11*v4
        dwdz = b3*w1 + b6*w2 + b9*w3 + b12*w4
        dudy = b2*u1 + b5*u2 + b8*u3 + b11*u4
        dvdx = b1*v1 + b4*v2 + b7*v3 + b10*v4
        dvdz = b3*v1 + b6*v2 + b9*v3 + b12*v4
        dwdy = b2*w1 + b5*w2 + b8*w3 + b11*w4
        dudz = b3*u1 + b6*u2 + b9*u3 + b12*u4
        dwdx = b1*w1 + b4*w2 + b7*w3 + b10*w4
! 
! deformation gradients F
!
        f11 = (1.d0 + dudx)
        f12 = dudy
        f13 = dudz
        f21 = dvdx
        f22 = (1.d0 + dvdy)
        f23 = dvdz
        f31 = dwdx
        f32 = dwdy
        f33 = (1.d0 + dwdz)

!             4
!            ----   ,
! Evaluate   \     P 
!            /      a
!            ----
!            a = 1

        s11(1) = 0.d0
        s22(1) = 0.d0
        s33(1) = 0.d0
        s12(1) = 0.d0
        s23(1) = 0.d0
        s13(1) = 0.d0
        
        DO k = 1, 4
           s11(1) = s11(1) + alpha(k,ielnum)*s11n(nodes(k,ielnum))
           s22(1) = s22(1) + alpha(k,ielnum)*s22n(nodes(k,ielnum))
           s33(1) = s33(1) + alpha(k,ielnum)*s33n(nodes(k,ielnum))
           s12(1) = s12(1) + alpha(k,ielnum)*s12n(nodes(k,ielnum))
           s23(1) = s23(1) + alpha(k,ielnum)*s23n(nodes(k,ielnum))
           s13(1) = s13(1) + alpha(k,ielnum)*s13n(nodes(k,ielnum)) 
        ENDDO

! calculate the volume

        vol = vx6 / 6.d0

!!$
!!$! ASSEMBLE THE INTERNAL FORCE VECTOR
!!$!
!!$! local node 1
!!$
!!$         Rnet(k1i) = Rnet(k1i) - VolEl0* &
!!$              ( S11e*B1*(1.d0+dudx) + S22e*B2*dudy + S33e*B3*dudz &
!!$              + S12e*( B2*(1.d0+dudx) + B1*dudy ) &
!!$              + S23e*( B3*dudy + B2*dudz ) &
!!$              + S13e*( B3*(1.d0+dudx) + B1*dudz ) )
!!$            Rnet(k2i) = Rnet(k2i) - VolEl0* &
!!$                 ( S11e*B1*dvdx + S22e*B2*(1.d0+dvdy) + S33e*B3*dvdz &
!!$                 + S12e*( B1*(1.d0+dvdy) + B2*dvdx ) &
!!$                 + S23e*( B3*(1.d0+dvdy) + B2*dvdz ) &
!!$                 + S13e*( B3*dvdx + B1*dvdz ) )
!!$            Rnet(k3i)   = Rnet(k3i)   - VolEl0* &
!!$                 ( S11e*B1*dwdx + S22e*B2*dwdy + S33e*B3*(1.d0+dwdz) &
!!$                 + S12e*( B2*dwdx + B1*dwdy ) &
!!$                 + S23e*( B3*dwdy + B2*(1.d0 + dwdz) ) &
!!$                 + S13e*( B3*dwdx + B1*(1.d0 + dwdz) ) ) 
!!$! local node 2 
!!$            Rnet(k1i) = Rnet(k1i) - VolEl0* &
!!$                 ( S11e*B4*(1.d0+dudx) + S22e*B5*dudy + S33e*B6*dudz &
!!$                 + S12e*( B5*(1.d0+dudx) + B4*dudy ) &
!!$                 + S23e*( B6*dudy + B5*dudz ) &
!!$                 + S13e*( B6*(1.d0+dudx) + B4*dudz ) ) 
!!$            Rnet(k2i) = Rnet(k2i) - VolEl0* &
!!$                 ( S11e*B4*dvdx + S22e*B5*(1.d0+dvdy) + S33e*B6*dvdz &
!!$                 + S12e*( B4*(1.d0+dvdy) + B5*dvdx ) &
!!$                 + S23e*( B6*(1.d0+dvdy) + B5*dvdz ) &
!!$                 + S13e*( B6*dvdx + B4*dvdz ) )
!!$            Rnet(k3i)   = Rnet(k3i)   - VolEl0* &
!!$                 ( S11e*B4*dwdx + S22e*B5*dwdy + S33e*B6*(1.d0+dwdz) &
!!$                 + S12e*( B5*dwdx + B4*dwdy ) &
!!$                 + S23e*( B6*dwdy + B5*(1.d0 + dwdz) ) &
!!$                 + S13e*( B6*dwdx + B4*(1.d0 + dwdz) ) )
!!$! local node 3 
!!$            Rnet(k1i) = Rnet(k1i) - VolEl0* &
!!$                 ( S11e*B7*(1.d0+dudx) + S22e*B8*dudy + S33e*B9*dudz &
!!$                 + S12e*( B8*(1.d0+dudx) + B7*dudy ) &
!!$                 + S23e*( B9*dudy + B8*dudz ) &
!!$                 + S13e*( B9*(1.d0+dudx) + B7*dudz ) )
!!$            Rnet(k2i) = Rnet(k2i) - VolEl0* &
!!$                 ( S11e*B7*dvdx + S22e*B8*(1.d0+dvdy) + S33e*B9*dvdz &
!!$                 + S12e*( B7*(1.d0+dvdy) + B8*dvdx ) &
!!$                 + S23e*( B9*(1.d0+dvdy) + B8*dvdz ) &
!!$                 + S13e*( B9*dvdx + B7*dvdz ) )
!!$            Rnet(k3i)   = Rnet(k3i)   - VolEl0* &
!!$                 ( S11e*B7*dwdx + S22e*B8*dwdy + S33e*B9*(1.d0+dwdz) &
!!$                 + S12e*( B8*dwdx + B7*dwdy ) &
!!$                 + S23e*( B9*dwdy + B8*(1.d0 + dwdz) ) &
!!$                 + S13e*( B9*dwdx + B7*(1.d0 + dwdz) ) )
!!$! local node 4         
!!$            Rnet(k1i) = Rnet(k1i) - VolEl0* &
!!$                 ( S11e*B10*(1.d0+dudx) + S22e*B11*dudy+S33e*B12*dudz &
!!$                 + S12e*( B11*(1.d0+dudx) + B10*dudy ) &
!!$                 + S23e*( B12*dudy + B11*dudz ) &
!!$                 + S13e*( B12*(1.d0+dudx) + B10*dudz ) )
!!$            Rnet(k2i) = Rnet(k2i) - VolEl0* &
!!$                 ( S11e*B10*dvdx + S22e*B11*(1.d0+dvdy)+S33e*B12*dvdz &
!!$                 + S12e*( B10*(1.d0+dvdy) + B11*dvdx ) &
!!$                 + S23e*( B12*(1.d0+dvdy) + B11*dvdz ) &
!!$                 + S13e*( B12*dvdx + B10*dvdz ) ) 
!!$            Rnet(k3i)   = Rnet(k3i)   - VolEl0* &
!!$                 ( S11e*B10*dwdx + S22e*B11*dwdy+S33e*B12*(1.d0+dwdz) &
!!$                 + S12e*( B11*dwdx + B10*dwdy ) &
!!$                 + S23e*( B12*dwdy + B11*(1.d0 + dwdz) ) &
!!$                 + S13e*( B12*dwdx + B10*(1.d0 + dwdz) ) )

        rnet(k1n1) = rnet(k1n1) - volel0* &
             ( s11(1)*b1*f11 + s22(1)*b2*f12 + s33(1)*b3*f13 &
             + s12(1)*( b2*f11 + b1*f12 ) &
             + s23(1)*( b3*f12 + b2*f13 ) &
             + s13(1)*( b3*f11 + b1*f13 ) )
        rnet(k2n1) = rnet(k2n1) - volel0* &
             ( s11(1)*b1*f21 + s22(1)*b2*f22 + s33(1)*b3*f23 &
             + s12(1)*( b1*f22 + b2*f21 ) &
             + s23(1)*( b3*f22 + b2*f23 ) &
             + s13(1)*( b3*f21 + b1*f23 ) )
        rnet(k3n1) = rnet(k3n1) - volel0* &
             ( s11(1)*b1*f31 + s22(1)*b2*f32 + s33(1)*b3*f33 &
             + s12(1)*( b2*f31 + b1*f32 ) &
             + s23(1)*( b3*f32 + b2*f33 ) &
             + s13(1)*( b3*f31 + b1*f33 ) )
! local node 2 
        rnet(k1n2) = rnet(k1n2) - volel0* &
             ( s11(1)*b4*f11 + s22(1)*b5*f12 + s33(1)*b6*f13 &
             + s12(1)*( b5*f11 + b4*f12 ) &
             + s23(1)*( b6*f12 + b5*f13 ) &
             + s13(1)*( b6*f11 + b4*f13 ) )
        rnet(k2n2) = rnet(k2n2) - volel0* &
             ( s11(1)*b4*f21 + s22(1)*b5*f22 + s33(1)*b6*f23 &
             + s12(1)*( b4*f22 + b5*f21 ) &
             + s23(1)*( b6*f22 + b5*f23 ) &
             + s13(1)*( b6*f21 + b4*f23 ) )
        rnet(k3n2) = rnet(k3n2) - volel0* &
             ( s11(1)*b4*f31 + s22(1)*b5*f32 + s33(1)*b6*f33 &
             + s12(1)*( b5*f31 + b4*f32 ) &
             + s23(1)*( b6*f32 + b5*f33 ) &
             + s13(1)*( b6*f31 + b4*f33 ) )
! local node 3 
        rnet(k1n3) = rnet(k1n3) - volel0* &
             ( s11(1)*b7*f11 + s22(1)*b8*f12 + s33(1)*b9*f13 &
             + s12(1)*( b8*f11 + b7*f12 ) &
             + s23(1)*( b9*f12 + b8*f13 ) &
             + s13(1)*( b9*f11 + b7*f13 ) )
        rnet(k2n3) = rnet(k2n3) - volel0* &
             ( s11(1)*b7*f21 + s22(1)*b8*f22 + s33(1)*b9*f23 &
             + s12(1)*( b7*f22 + b8*f21 ) &
             + s23(1)*( b9*f22 + b8*f23 ) &
             + s13(1)*( b9*f21 + b7*f23 ) )
        rnet(k3n3) = rnet(k3n3) - volel0* &
             ( s11(1)*b7*f31 + s22(1)*b8*f32 + s33(1)*b9*f33 &
             + s12(1)*( b8*f31 + b7*f32 ) &
             + s23(1)*( b9*f32 + b8*f33 ) &
             + s13(1)*( b9*f31 + b7*f33 ) )
! local node 4         
        rnet(k1n4) = rnet(k1n4) - volel0* &
             ( s11(1)*b10*f11 + s22(1)*b11*f12+s33(1)*b12*f13 &
             + s12(1)*( b11*f11 + b10*f12 ) &
             + s23(1)*( b12*f12 + b11*f13 ) &
             + s13(1)*( b12*f11 + b10*f13 ) )
        rnet(k2n4) = rnet(k2n4) - volel0* &
             ( s11(1)*b10*f21 + s22(1)*b11*f22 + s33(1)*b12*f23 &
             + s12(1)*( b10*f22 + b11*f21 ) &
             + s23(1)*( b12*f22 + b11*f23 ) &
             + s13(1)*( b12*f21 + b10*f23 ) )
        rnet(k3n4) = rnet(k3n4) - volel0* &
             ( s11(1)*b10*f31 + s22(1)*b11*f32 + s33(1)*b12*f33 &
             + s12(1)*( b11*f31 + b10*f32 ) &
             + s23(1)*( b12*f32 + b11*f33 ) &
             + s13(1)*( b12*f31 + b10*f33 ) )
        
     ENDDO
  ENDDO

  RETURN
END SUBROUTINE v3d4n_neohookeancompress