v3d4n_nl.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(numnp, numel, coor, disp, nodes, Rnet,&
     numelneigh, elconn, alpha, ahat,nummatvol, &
     xmu,xlambda)

!!****f* Rocfrac/Source/v3d4n_nl.f90
!!
!!  NAME
!!    v3d4n_nl
!!
!!  FUNCTION
!!
!!    Computes the internal force vector for a 4-node tetrahedral
!!    NODE BASED ELEMENT  for finite deformations.
!!
!!  INPUTS
!!    
!!    NumNP -- Number of nodes
!!   NumEL  -- Number of elements
!!    coor  -- number of coordinates
!!   disp   -- Nodal Displacement
!!   nodes --  Nodal connectivity
!!   NumElNeigh -- Number of elements in contact with node
!!   ElConn -- Element connectivity
!!   alpha -- volume ratio
!!   Ahat -- undeformed volume of node
!!   NumMatVol -- number of materials
!!   xmu -- material parameter
!!   xlambda -- material parameter
!!  
!!  OUTPUT
!!
!!    Rnet -- internal force vector
!!
!!****

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

! Fix, not parallel


  INTEGER :: i,k,j
  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
!--   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
  REAL*8 :: x12, x13, y12, y13
!-- Added these to speed up B calculation

  REAL*8 ::  z12, z13
  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 :: 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 :: vx6def, vx6invdef, voldef, vtotal

  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
     
     f11 = f11*vainv
     f22 = f22*vainv
     f33 = f33*vainv
     f12 = f12*vainv
     f13 = f13*vainv
     f21 = f21*vainv
     f23 = f23*vainv
     f31 = f31*vainv
     f32 = f32*vainv
! should this be before divide by volume
     jacob = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)

!       Jacob = voldef*VaInv
     
     IF( jacob.LE.0.d0) THEN
        WRITE(*,100) i
        stop
     ENDIF
!
! Compressible Neo-Hookean Material
! Eq. (5.29)
       
     coeff1 = xmu(1)/jacob
     coeff2 = xlambda(1)/jacob*log(jacob)
     
     s11n(i) = coeff1*(f11**2+f12**2+f13**2-1.d0) + coeff2
     s12n(i) = coeff1*(f21*f11+f12*f22+f13*f23)
     s13n(i) = coeff1*(f31*f11+f12*f32+f13*f33)
     s22n(i) = coeff1*(f21**2+f22**2+f23**2-1.d0) + coeff2
     s23n(i) = coeff1*(f21*f31+f32*f22+f23*f33)
     s33n(i) = coeff1*(f31**2+f32**2+f33**2-1.d0) + coeff2
     
  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
     
     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