MP/Source/v3d4n_nl_arrudaboyce.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_nl_arruda_boyce(numnp, numel, coor, disp, nodes, Rnet, &
     s11ev, s22ev, s33ev, s12ev, s23ev, s13ev, &
     numelneigh,elconn,alpha, ahat,&
     xmu,xkappa,nummatvol, &
     nprocs,totnumndcomm,totnumneighprocs,neighproclist,numndcomm,neigh_lst,&
     mpi_status_size,mpi_comm_rocfrac,mpi_double_precision, &
     reqrcv,reqsnd,statrcv, statsnd)

    USE rocstar_rocfraccomm 

    IMPLICIT NONE

    INTEGER :: k1, j1, k2

    INTEGER :: totnumneighprocs, nprocs,totnumndcomm
    INTEGER, DIMENSION(1:TotNumNeighProcs) :: neighproclist
    INTEGER, DIMENSION(1:TotNumNeighProcs) ::numndcomm

    TYPE(rcv_buf), pointer, DIMENSION(:) :: recvdatafrm 
    TYPE(send_buf), DIMENSION(1:TotNumNeighProcs) :: neigh_lst  

    integer :: nummatvol
    INTEGER :: numnp, numel
    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:NumMatVol) :: xmu, xkappa
    REAL*8, DIMENSION(1:numel) :: s11ev, s22ev, s33ev, s12ev, s23ev, s13ev
    
    INTEGER :: j,i,k
    REAL*8 :: aaa
!--   coordinate holding variable
    REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
    REAL*8 :: x1d,x2d,x3d,x4d,y1d,y2d,y3d,y4d,z1d,z2d,z3d,z4d
    INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
    INTEGER :: k3n1,k3n2,k3n3,k3n4

    REAL*8 :: ahatinv
    INTEGER :: ielnum
    REAL*8 :: sixinv
!--   x, y, and z displacements of nodes
    REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
!--   partial derivatives of the displacement 
    REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx

!--  Coordinate subtractions
    REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!-- Added these to speed up B calculation
    
    REAL*8 ::  z12, z13,x12,x13,y12,y13
    REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33
    REAL*8 :: sh1, sh2, sh3, sh4, sh5, sh6, sh7, sh8, sh9, sh10, sh11, sh12
    REAL*8 :: sh13, sh14, sh15
!   --   6*volume and the volume      
    REAL*8 :: vx6
!--   spacial derivatives
    REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
!--   strains
    REAL*8 :: e11,e22,e33,e12,e23,e13
    REAL*8 :: b13, b14, b15
    REAL*8 :: sums11n,sums22n,sums33n,sums12n,sums23n,sums13n
    REAL*8 :: sume11n,sume22n,sume33n,sume12n,sume23n,sume13n
    REAL*8, DIMENSION(1:numnp) :: s11n, s22n, s33n,s12n, s23n, s13n
    
    REAL*8 :: e11n, e22n, e33n, e12n, e23n, e13n
    INTEGER :: k3i,k2i,k1i
    INTEGER :: ix
    INTEGER :: n1, n2, n3, n4
    REAL*8 :: volel, volel0,vx6inv,vol
    REAL*8 :: s11e, s22e, s33e, s12e, s23e, s13e
    REAL*8 :: vainv

    REAL*8 :: coeff1, coeff2
    REAL*8 :: jacob
    REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
    REAL*8,DIMENSION(1:numnp) :: f11v, f12v, f13v, f21v, f22v, f23v, f31v, f32v, f33v, voldefv
    REAL*8 :: vx6def, vx6invdef, voldef


    INTEGER :: myid
    INTEGER :: ierr

    REAL*8, pointer, dimension(:) :: buf

    integer :: mpi_status_size,mpi_comm_rocfrac,mpi_double_precision

!--   Non-block receive, Non-block send request arrays
    INTEGER, DIMENSION(1:TotNumNeighProcs) :: reqrcv, reqsnd
    INTEGER, DIMENSION(1:MPI_STATUS_SIZE,1:TotNumNeighProcs) :: statrcv, statsnd

    sixinv = 1.d0/6.d0

    DO i = 1, numnp ! for each node

! /* Undeformed Volume of node */

       vainv = 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 

       voldef = 0.d0
              
       DO j = 1, numelneigh(i) ! for each element associated with node 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) 

! 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
! 
! deformation gradients F
!
          f11 = f11 + alpha(ix,ielnum)*volel0*(1.d0 + ( b1*u1 + b4*u2 + b7*u3 + b10*u4 )) ! 1 + ( dudx )
          f22 = f22 + alpha(ix,ielnum)*volel0*(1.d0 + ( b2*v1 + b5*v2 + b8*v3 + b11*v4 )) ! 1 + ( dvdy )
          f33 = f33 + alpha(ix,ielnum)*volel0*(1.d0 + ( b3*w1 + b6*w2 + b9*w3 + b12*w4 )) ! 1 + ( dwdz )
          f12 = f12 + alpha(ix,ielnum)*volel0*(b2*u1 + b5*u2 + b8*u3 + b11*u4) ! dudy
          f21 = f21 + alpha(ix,ielnum)*volel0*(b1*v1 + b4*v2 + b7*v3 + b10*v4) ! dvdx
          f23 = f23 + alpha(ix,ielnum)*volel0*(b3*v1 + b6*v2 + b9*v3 + b12*v4) ! dvdz
          f32 = f32 + alpha(ix,ielnum)*volel0*(b2*w1 + b5*w2 + b8*w3 + b11*w4) ! dwdy
          f13 = f13 + alpha(ix,ielnum)*volel0*(b3*u1 + b6*u2 + b9*u3 + b12*u4) ! dudz
          f31 = f31 + alpha(ix,ielnum)*volel0*(b1*w1 + b4*w2 + b7*w3 + b10*w4) ! dwdx

       ENDDO

       f11v(i) = f11*vainv
       f22v(i) = f22*vainv
       f33v(i) = f33*vainv
       f12v(i) = f12*vainv
       f13v(i) = f13*vainv
       f21v(i) = f21*vainv
       f23v(i) = f23*vainv
       f31v(i) = f31*vainv
       f32v(i) = f32*vainv

       voldefv(i) = voldef

    ENDDO



    ALLOCATE(recvdatafrm(0:nprocs-1))
!
!----- FORM THE BUFFER CONTAINING COMMUNICATED STRESS MATRIX NODAL VALUES
    
    ALLOCATE(buf(1:totnumndcomm/3*10)) ! is really (TotNumNdComm*3 * 3 + 1)
    k1 = 1
    DO j1 = 1, totnumneighprocs
       k = neighproclist(j1)
       ALLOCATE(recvdatafrm(k)%rcvbuf(1:numndcomm(j1)*10))
       DO j = 1, numndcomm(j1)
          k2 = neigh_lst(j1)%NdId(j) !NdCommList(k)%NdId(j)
          buf(k1  ) = f11v(k2)
          buf(k1+1) = f12v(k2)
          buf(k1+2) = f13v(k2)
          buf(k1+3) = f21v(k2)
          buf(k1+4) = f22v(k2)
          buf(k1+5) = f23v(k2)
          buf(k1+6) = f31v(k2)
          buf(k1+7) = f32v(k2)
          buf(k1+8) = f33v(k2)
          buf(k1+9) = voldefv(k2)
          k1 = k1 + 10
       ENDDO
    ENDDO
!
!-MPI- RECEIVE THE RECIPRICAL MASS MATRIX DIAGONAL FROM THE NEIGHBORS
!
    DO j1 = 1, totnumneighprocs
       k = neighproclist(j1)
       CALL mpi_irecv(recvdatafrm(k)%rcvbuf(1),numndcomm(j1)*10, &
            mpi_double_precision, k, 10, mpi_comm_rocfrac,reqrcv(j1),ierr)
    ENDDO
!
!-MPI- SEND THE RECIPRICAL MASS MATRIX DIAGONAL TO THE NEIGHBORS
!
    k2 = 1
    DO j1 = 1, totnumneighprocs
       k = neighproclist(j1)
       
       CALL mpi_isend(buf(k2),numndcomm(j1)*10,&
            mpi_double_precision,k,10,mpi_comm_rocfrac,reqsnd(j1),ierr)
       k2 = k2 + numndcomm(j1)*10
    ENDDO
!
!-MPI- WAIT FOR INTERNAL FORCE VECTOR COMMUNICATION TO COMPLETE
!
    IF(totnumneighprocs.GT.0)THEN
       CALL mpi_waitall(totnumneighprocs,reqrcv,statrcv,ierr)
       CALL mpi_waitall(totnumneighprocs,reqsnd,statsnd,ierr)
    ENDIF
    DEALLOCATE(buf)

!
!----- ADD NEIGHBOR'S CONTRIBUTION TO THE INTERNAL FORCE VECTOR
!
    DO j1 = 1, totnumneighprocs
       k = neighproclist(j1)
       k1 = 1
       DO j = 1, numndcomm(j1)
          k2 = neigh_lst(j1)%NdId(j)
          f11v(k2)  = f11v(k2)  + recvdatafrm(k)%rcvbuf(k1)
          f12v(k2)  = f12v(k2)  + recvdatafrm(k)%rcvbuf(k1+1)
          f13v(k2)  = f13v(k2)  + recvdatafrm(k)%rcvbuf(k1+2)
          f21v(k2)  = f21v(k2)  + recvdatafrm(k)%rcvbuf(k1+3)
          f22v(k2)  = f22v(k2)  + recvdatafrm(k)%rcvbuf(k1+4)
          f23v(k2)  = f23v(k2)  + recvdatafrm(k)%rcvbuf(k1+5)
          f31v(k2)  = f31v(k2)  + recvdatafrm(k)%rcvbuf(k1+6)
          f32v(k2)  = f32v(k2)  + recvdatafrm(k)%rcvbuf(k1+7)
          f33v(k2)  = f33v(k2)  + recvdatafrm(k)%rcvbuf(k1+8)
          voldefv(k2) = voldefv(k2)  + recvdatafrm(k)%rcvbuf(k1+9)
          k1 = k1 + 10
       ENDDO
    ENDDO

    DEALLOCATE(recvdatafrm)
    DO i = 1, numnp ! for each node

       f11 = f11v(i)
       f22 = f22v(i)
       f33 = f33v(i)
       f12 = f12v(i)
       f13 = f13v(i)
       f21 = f21v(i)
       f23 = f23v(i)
       f31 = f31v(i)
       f32 = f32v(i)

       voldef = voldefv(i)

! WRONG FOR MORE THEN ONE MATERIAL

       CALL arruda_boyce_cauchy(f11,f12,f13,f21,f22,f23,f31,f32,f33, &
            s11n(i),s22n(i),s33n(i),s12n(i),s13n(i),s23n(i),i,xmu(1),xkappa(1))
    ENDDO


    DO ielnum = 1, numel  ! For each element

       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)
       
       x1 = x1 + u1
       x2 = x2 + u2
       x3 = x3 + u3
       x4 = x4 + u4
       y1 = y1 + v1
       y2 = y2 + v2
       y3 = y3 + v3
       y4 = y4 + v4
       z1 = z1 + w1
       z2 = z2 + w2
       z3 = z3 + w3
       z4 = z4 + w4

       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
       
       vx6def = -( x14*c11 + y14*c21 + z14*c31 )
       
       vx6invdef = 1.d0 / vx6def 

! calculate the volume
       voldef = vx6def/6.d0
            
       IF(vx6def.LE.0.d0) THEN
          WRITE(*,100) i
          stop
       ENDIF

    ! equation (14)

       s11e = 0.d0
       s22e = 0.d0
       s33e = 0.d0
       s12e = 0.d0
       s23e = 0.d0
       s13e = 0.d0
       
       DO k = 1, 4
          s11e = s11e + alpha(k,ielnum)*s11n(nodes(k,ielnum))
          s22e = s22e + alpha(k,ielnum)*s22n(nodes(k,ielnum))
          s33e = s33e + alpha(k,ielnum)*s33n(nodes(k,ielnum))
          s12e = s12e + alpha(k,ielnum)*s12n(nodes(k,ielnum))
          s23e = s23e + alpha(k,ielnum)*s23n(nodes(k,ielnum))
          s13e = s13e + alpha(k,ielnum)*s13n(nodes(k,ielnum))   
       ENDDO

       s11ev(ielnum) = s11e
       s22ev(ielnum) = s22e
       s33ev(ielnum) = s33e
       s12ev(ielnum) = s12e
       s23ev(ielnum) = s23e
       s13ev(ielnum) = s13e

       sh1  = (y34*z24 - y24*z34) * vx6invdef
       sh2  = (z34*x24 - z24*x34) * vx6invdef
       sh3  = (x34*y24 - x24*y34) * vx6invdef
       sh4  = (y13*z14 - y14*z13) * vx6invdef
       sh5  = (z13*x14 - z14*x13) * vx6invdef
       sh6  = (x13*y14 - x14*y13) * vx6invdef
       sh7  = (y14*z12 - y12*z14) * vx6invdef
       sh8  = (z14*x12 - z12*x14) * vx6invdef
       sh9  = (x14*y12 - x12*y14) * vx6invdef
       sh10 = (y12*z13 - y13*z12) * vx6invdef
       sh11 = (z12*x13 - z13*x12) * vx6invdef
       sh12 = (x12*y13 - x13*y12) * vx6invdef

! ASSEMBLE THE INTERNAL FORCE VECTOR
!
! local node 1

       rnet(k1n1) = rnet(k1n1) - voldef* &
            ( s11e*sh1 + s12e*sh2 + s13e*sh3 )
       rnet(k2n1) = rnet(k2n1) - voldef* &
            ( s12e*sh1 + s22e*sh2 + s23e*sh3 )
       rnet(k3n1) = rnet(k3n1) - voldef* &
            ( s13e*sh1 + s23e*sh2 + s33e*sh3 )
! local node 2 
       rnet(k1n2) = rnet(k1n2) - voldef* &
            ( s11e*sh4 + s12e*sh5 + s13e*sh6 )
       rnet(k2n2) = rnet(k2n2) - voldef* &
            ( s12e*sh4 + s22e*sh5 + s23e*sh6 )
       rnet(k3n2) = rnet(k3n2) - voldef* &
            ( s13e*sh4 + s23e*sh5 + s33e*sh6 )
! local node 3 
       rnet(k1n3) = rnet(k1n3) - voldef* &
            ( s11e*sh7 + s12e*sh8 + s13e*sh9 )
       rnet(k2n3) = rnet(k2n3) - voldef* &
            ( s12e*sh7 + s22e*sh8 + s23e*sh9 )
       rnet(k3n3) = rnet(k3n3) - voldef* &
            ( s13e*sh7 + s23e*sh8 + s33e*sh9 )
! local node 4  
       rnet(k1n4) = rnet(k1n4) - voldef* &
            ( s11e*sh10 + s12e*sh11 + s13e*sh12 )
       rnet(k2n4) = rnet(k2n4) - voldef* &
            ( s12e*sh10 + s22e*sh11 + s23e*sh12 )
       rnet(k3n4) = rnet(k3n4) - voldef* &
            ( s13e*sh10 + s23e*sh11 + s33e*sh12 )
    
    ENDDO

    RETURN
100 FORMAT(' Negative Jacobian for element: ',i10)

  END SUBROUTINE v3d4n_nl_arruda_boyce