TractPressLoad.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 tractpressload(R_ex,numnp, &
     interfacesfnumelems, interfacesfnumnodes, &
     interfacesfelemconn, &
     mapnodesf,lwrbnd,uppbnd, coor, disp, mapsfelvolel,&
     elconnvol,ieltype,numelvol,tractpress )

  IMPLICIT NONE

! Transforms the pressure given in the deformed state to
! the undeformed configuration. This subroutine is for
! the formulation where all quantities are with respect to
! the undeformed configuration.  Assumes the 
! pressure (tractions) are constant over the surface of
! the triangle.

! 1-3 for 4 node tet (i.e. 3 node triangles)
! 4-6 for 10 node tet (i.e. 6 node triangles, mid-side nodes get traction)

  INTEGER :: lwrbnd,uppbnd

  INTEGER :: ieltype,numelvol
  INTEGER, DIMENSION(1:iElType,1:NumElVol) :: elconnvol
  
  INTEGER :: numnp          ! number of nodes

  INTEGER :: interfacesfnumelems,interfacesfnumnodes
  INTEGER, DIMENSION(1:UppBnd,1:InterfaceSFNumElems) :: interfacesfelemconn
  INTEGER, DIMENSION(1:InterfaceSFNumNodes) :: mapnodesf
  REAL*8, DIMENSION(1:3*NumNp) :: disp
  
  REAL*8, DIMENSION(1:3*numnp) :: r_ex ! external force
  
  INTEGER :: nx, ny, nz
  INTEGER :: i,j,k
  
  INTEGER,DIMENSION(1:3) :: triconn
  REAL*8,DIMENSION(1:3)  :: uniftract

! mesh coordinates
  REAL*8, DIMENSION(1:3,1:numnp) :: coor

  REAL*8 :: x0p,y0p,z0p
  REAL*8 :: x1p,y1p,z1p
  REAL*8 :: x2p,y2p,z2p
  REAL*8 :: x3p,y3p,z3p

! Components sides vector

  REAL*8 :: x1x0, y1y0, z1z0
  REAL*8 :: x2x0, y2y0, z2z0
  
  REAL*8 :: magnorm, area, at
  REAL*8 :: gauss =  0.333333333333333d0
  REAL*8 :: weight = 1.0d0
  REAL*8 :: xnorm,ynorm,znorm
  REAL*8 :: u1,v1,w1,u2,v2,w2,u3,v3,w3,u4,v4,w4
  REAL*8 :: x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4
  REAL*8 :: j11,j12,j13,j21,j22,j23
  REAL*8 :: invdenjplus
  REAL*8 :: jplus11, jplus12, jplus21, jplus22, jplus31, jplus32
  REAL*8 :: gradu11, gradu12, gradu13, gradu21, gradu22, gradu23
  REAL*8 :: gradu31, gradu32, gradu33
  REAL*8 :: f11t,f12t,f13t,f21t,f22t,f23t,f31t,f32t,f33t
  REAL*8 :: j1, undefnorm(1:3)
  INTEGER :: idvol

  INTEGER :: n1,n2,n3,n4
! -- dummy and counters
  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
  INTEGER :: k3n1,k3n2,k3n3,k3n4
  REAL*8 :: tractpress
!--  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
  REAL*8 :: vx6, vx6inv
!--   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
  INTEGER, DIMENSION(1:NumElVol) :: mapsfelvolel
  REAL*8 :: jac3d
  REAL*8 :: defnorm(1:3)
  REAL*8 :: denfinvt,magnormdef,area0
  REAL*8 :: f11,f12,f13,f21,f22,f23,f31,f32,f33
  REAL*8 :: finvt11, finvt12, finvt13, finvt21, finvt22, finvt23, finvt31, finvt32, finvt33

  INTEGER :: n1s,n2s,n3s

  DO i = 1, interfacesfnumelems

! Compute the Gradient of U (disp) Volumetric 

     idvol = mapsfelvolel(i) ! Relates Surface id to Volumetric id
     
     n1 = elconnvol(1,idvol)
     n2 = elconnvol(2,idvol)
     n3 = elconnvol(3,idvol)
     n4 = elconnvol(4,idvol)

     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

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

     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

! See the maple worksheet 'V3D4.mws' for the derivation of [B]
! 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

!
!  
!

     denfinvt = f11*(f22*f33-f32*f23) - f12*(f21*f33+f31*f23) + f13*(f21*f32-f31*f22)

     finvt11 =  (f22*f33-f32*f23)/denfinvt
     finvt12 = -(f21*f33-f31*f23)/denfinvt
     finvt13 =  (f21*f32-f31*f22)/denfinvt
     finvt21 = -(f12*f33-f32*f13)/denfinvt
     finvt22 =  (f11*f33-f31*f13)/denfinvt
     finvt23 = -(f11*f32-f31*f12)/denfinvt
     finvt31 =  (f12*f23-f22*f13)/denfinvt
     finvt32 = -(f11*f23-f21*f13)/denfinvt
     finvt33 =  (f11*f22-f21*f12)/denfinvt

     jac3d = denfinvt 
     
!  Nodes of the Triangle ( Vertex Nodes ) 

     triconn(1:3) = interfacesfelemconn(1:3,i)

! 1) Calculate the Normal in the Deformed Configuration

!  Uses the fact that the norm of the cross product vector
!  is the area of the parallelogram they form.  The triangle they
!  form has half that area.

     n1s = mapnodesf(triconn(1))
     n2s = mapnodesf(triconn(2))
     n3s = mapnodesf(triconn(3))

     !TEST
     IF(n1s.NE.n1.AND.n1s.NE.n2.AND.n1s.NE.n3.AND.n1s.NE.n4)THEN
        print*,'n1s NOT part of volume element',n1s,n1,n2,n3,n4,mapsfelvolel(i),i
        stop
     ENDIF
     IF(n2s.NE.n1.AND.n2s.NE.n2.AND.n2s.NE.n3.AND.n2s.NE.n4)THEN
        print*,'n2s NOT part of volume element'
        stop
     ENDIF
     IF(n3s.NE.n1.AND.n3s.NE.n2.AND.n3s.NE.n3.AND.n3s.NE.n4)THEN
        print*,'n3s NOT part of volume element'
        stop
     ENDIF
        

     x1 = coor(1,mapnodesf(triconn(1)))
     y1 = coor(2,mapnodesf(triconn(1)))
     z1 = coor(3,mapnodesf(triconn(1)))
     
     x2 = coor(1,mapnodesf(triconn(2)))
     y2 = coor(2,mapnodesf(triconn(2)))
     z2 = coor(3,mapnodesf(triconn(2)))
     
     x3 = coor(1,mapnodesf(triconn(3)))
     y3 = coor(2,mapnodesf(triconn(3)))
     z3 = coor(3,mapnodesf(triconn(3)))

     u1 = disp(3*mapnodesf(triconn(1))-2)
     v1 = disp(3*mapnodesf(triconn(1))-1)
     w1 = disp(3*mapnodesf(triconn(1))  )
     
     u2 = disp(3*mapnodesf(triconn(2))-2)
     v2 = disp(3*mapnodesf(triconn(2))-1)
     w2 = disp(3*mapnodesf(triconn(2))  )
     
     u3 = disp(3*mapnodesf(triconn(3))-2)
     v3 = disp(3*mapnodesf(triconn(3))-1)
     w3 = disp(3*mapnodesf(triconn(3))  )

     x0p = x1
     y0p = y1
     z0p = z1 
     
     x1p = x2
     y1p = y2 
     z1p = z2
     
     x2p = x3
     y2p = y3
     z2p = z3


! Find vector componets of the sides

     x1x0 = x1p - x0p
     y1y0 = y1p - y0p
     z1z0 = z1p - z0p

     x2x0 = x2p - x0p
     y2y0 = y2p - y0p
     z2z0 = z2p - z0p     

!  Take the cross product

     x3p = y1y0 * z2z0  -  z1z0  *  y2y0
     y3p = z1z0 * x2x0  -  x1x0  *  z2z0
     z3p = x1x0 * y2y0  -  y1y0  *  x2x0
         
!  Computes the Euclidean norm of a vector in 3D Deformed Config.
         
     magnorm = sqrt( x3p*x3p + y3p*y3p + z3p*z3p )

     xnorm = x3p/magnorm
     ynorm = y3p/magnorm
     znorm = z3p/magnorm

!     area0 = 0.5*MagNorm

! little n (deformed normal)

     defnorm(1) = jac3d * (finvt11*xnorm+finvt12*ynorm+finvt13*znorm)
     defnorm(2) = jac3d * (finvt21*xnorm+finvt22*ynorm+finvt23*znorm)
     defnorm(3) = jac3d * (finvt31*xnorm+finvt32*ynorm+finvt33*znorm)

     magnormdef = sqrt(  defnorm(1)**2 +  defnorm(2)**2 +  defnorm(3)**2 )

     defnorm(1) = defnorm(1)/magnormdef
     defnorm(2) = defnorm(2)/magnormdef
     defnorm(3) = defnorm(3)/magnormdef

     x0p = x1 + u1
     y0p = y1 + v1
     z0p = z1 + w1
     
     x1p = x2 + u2
     y1p = y2 + v2
     z1p = z2 + w2
     
     x2p = x3 + u3
     y2p = y3 + v3
     z2p = z3 + w3

! Find vector componets of the sides

     x1x0 = x1p - x0p
     y1y0 = y1p - y0p
     z1z0 = z1p - z0p

     x2x0 = x2p - x0p
     y2y0 = y2p - y0p
     z2z0 = z2p - z0p     

!  Take the cross product

     x3p = y1y0 * z2z0  -  z1z0  *  y2y0
     y3p = z1z0 * x2x0  -  x1x0  *  z2z0
     z3p = x1x0 * y2y0  -  y1y0  *  x2x0
         
!  Computes the Euclidean norm of a vector in 3D Deformed Config.

     magnormdef = sqrt( x3p*x3p + y3p*y3p + z3p*z3p )

     area = 0.5*magnormdef

! Pressure * Jac2D * n * Area0 (Jac2D = Area/Area0)

     uniftract(1) = -tractpress*defnorm(1)*area/3.d0
     uniftract(2) = -tractpress*defnorm(2)*area/3.d0
     uniftract(3) = -tractpress*defnorm(3)*area/3.d0

!     print*, -TractPress*defNorm(1)*Area/3.d0

!  Assembly into global external load vector R_ex

     triconn(1:3) = interfacesfelemconn(lwrbnd:uppbnd,i)

     DO j = 1, 3
        nz = mapnodesf(triconn(j))*3 ! volume node
        nx = nz - 2
        ny = nz - 1
        r_ex(nx) = r_ex(nx) + uniftract(1)
        r_ex(ny) = r_ex(ny) + uniftract(2)
        r_ex(nz) = r_ex(nz) + uniftract(3)
     ENDDO

  ENDDO

END SUBROUTINE tractpressload