TractLoad_Hex.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 tractload_hex(R_ex,numnp, &
     interfaceelemtract, &
     interfacenumelems, interfacenumnodes, &
     interfaceelemconn, &
     mapnode,lwrbnd,uppbnd,coor)

  IMPLICIT NONE

  INTEGER :: lwrbnd,uppbnd

  INTEGER :: interfacenumelems,interfacenumnodes
  INTEGER, DIMENSION(1:UppBnd,1:InterfaceNumElems) :: interfaceelemconn
  INTEGER, DIMENSION(1:InterfaceNumNodes) :: mapnode
  REAL*8,  DIMENSION(1:InterfaceNumElems) :: interfaceelemtract
  
  INTEGER :: numnp          ! number of nodes
  
  REAL*8, DIMENSION(1:3*numnp) :: r_ex ! external force
  
  INTEGER :: nx, ny, nz
  INTEGER :: i,j,k
  
  INTEGER,DIMENSION(1:4) :: conn2d
  real*8,DIMENSION(1:3) :: dlvec
  real*8 :: pmag, magnormdef, area

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

  

  REAL*8, DIMENSION(1:3) :: veca, vecb, outwd_surface_normal, &
       half_normal_veca, half_normal_vecb
  

  DO i = 1, interfacenumelems
     
!  Nodes of the Quad ( Vertex Nodes ) 

     conn2d(1:4) = interfaceelemconn(1:4,i)

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

     veca(1:3) =  coor(1:3,mapnode(conn2d(2)))- coor(1:3,mapnode(conn2d(1)))
     vecb(1:3) =  coor(1:3,mapnode(conn2d(4)))- coor(1:3,mapnode(conn2d(1)))
        
     half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
     half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
     half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
        
     veca(1:3) = coor(1:3,mapnode(conn2d(4))) - coor(1:3,mapnode(conn2d(3)))
     vecb(1:3) = coor(1:3,mapnode(conn2d(2))) - coor(1:3,mapnode(conn2d(3)))
        
     half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
     half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
     half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
        
     outwd_surface_normal = half_normal_veca + half_normal_vecb 

!     MagNormDef = SQRT( outwd_surface_normal(1)**2 + outwd_surface_normal(2)**2 + outwd_surface_normal(3)**2 )

!     Area = 0.5*MagNormDef

     pmag = interfaceelemtract(i)*0.125d0  ! /4 * /2
        
     dlvec(1:3) = - pmag * outwd_surface_normal(1:3)

! exactly the same results
!     print*,'**', - InterfaceElemTract(i) * Area * 0.25 * outwd_surface_normal(1:3)/MagNormDef
!     print*,'***', dlvec(1:3)


!  Assembly into global external load vector R_ex

     conn2d(1:4) = interfaceelemconn(lwrbnd:uppbnd,i)

     DO j = 1, 4
        nz = mapnode(conn2d(j))*3 ! volume node
        nx = nz - 2
        ny = nz - 1
        r_ex(nx) = r_ex(nx) + dlvec(1)
        r_ex(ny) = r_ex(ny) + dlvec(2)
        r_ex(nz) = r_ex(nz) + dlvec(3)
     ENDDO
  ENDDO
  
  RETURN
END SUBROUTINE tractload_hex