press_seg3.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    press_seg3            ../engine/source/loads/general/load_pcyl/press_seg3.F
!||--- called by ------------------------------------------------------
!||    h3d_pre_skin_scalar   ../engine/source/output/h3d/h3d_results/h3d_skin_scalar.F
!||    h3d_skin_vector       ../engine/source/output/h3d/h3d_results/h3d_skin_vector.F
!||    pressure_cyl          ../engine/source/loads/general/load_pcyl/pressure_cyl.F
!||--- calls      -----------------------------------------------------
!||    table_interp_dydx     ../engine/source/tools/curve/table_tools.F
!||--- uses       -----------------------------------------------------
!||    pload_cyl_mod         ../common_source/modules/loads/pload_cyl_mod.F
!||    table_mod             ../engine/share/modules/table_mod.F
!||====================================================================
      SUBROUTINE press_seg3(A     ,B     ,C     ,N1    ,DIR    ,
     .                      IFUNC ,TABLE ,XFACR ,XFACT ,PRESS  )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE table_mod
      USE pload_cyl_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER ,INTENT(IN) :: IFUNC
      my_real ,INTENT(IN) :: xfacr   !< radius scale factor
      my_real ,INTENT(IN) :: xfact   !< time   scale factor
      my_real ,DIMENSION(3) ,INTENT(IN)  :: a,b,c,n1,dir
      TYPE (TTABLE), DIMENSION(NTABLE) ,INTENT(IN)       :: TABLE
      my_real ,INTENT(OUT) :: press
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: NPT
      my_real :: a11,a12,a21,a22,b1,b2,det,alpha,beta,gamma,f1,f2,r,s,
     .   r1,r2,r3,rmax,p1,p2,p3,dydx,nx,ny,nz,norm,cosp,sinp
      my_real, DIMENSION(3) :: ab,ac,an,p,ap,bp,cp
      my_real, DIMENSION(2) :: xx
 
 
      INTEGER :: i
c-----------------------------------------------
c     1) find intersection point P between plane surface A-B-C and line N1-N2
c     2) calculate radial coordinates of each segment node in cylindrical
c        coord system defined by direction axis
c        (distance between point P and each segment node projected
c        on plane perpendicular to axis)
c=======================================================================
      ab(1) = b(1)  - a(1)
      ab(2) = b(2)  - a(2)
      ab(3) = b(3)  - a(3)
      ac(1) = c(1)  - a(1)
      ac(2) = c(2)  - a(2)
      ac(3) = c(3)  - a(3)
      an(1) = n1(1) - a(1)
      an(2) = n1(2) - a(2)
      an(3) = n1(3) - a(3)
      alpha = dir(1)
      beta  = dir(2)
      gamma = dir(3)
      npt   = SIZE(table(ifunc)%X(1)%VALUES)
      rmax  = table(ifunc)%X(1)%VALUES(npt) * xfacr
      p1    = zero
      p2    = zero
      p3    = zero
c
      IF (alpha /= zero) THEN
        f1 = beta  / alpha
        f2 = gamma / alpha
        b1  = an(2) - f1 * an(1)
        b2  = an(3) - f2 * an(1)
        a11 = ab(2) - f1 * ab(1)
        a12 = ac(2) - f1 * ac(1)
        a21 = ab(3) - f2 * ab(1)
        a22 = ac(3) - f2 * ac(1)
      ELSE IF (beta /= zero) THEN
        f1  = alpha / beta
        f2  = gamma / beta
        b1  = an(1) - f1 * an(2)
        b2  = an(3) - f2 * an(2)
        a11 = ab(1) - f1 * ab(2)
        a12 = ac(1) - f1 * ac(2)
        a21 = ab(3) - f2 * ab(2)
        a22 = ac(3) - f2 * ac(2)
      ELSE IF (gamma /= zero) THEN
        f1  = alpha / gamma
        f2  = beta  / gamma
        b1  = an(1) - f1 * an(3)
        b2  = an(2) - f2 * an(3)
        a11 = ab(1) - f1 * ab(3)
        a12 = ac(1) - f1 * ac(3)
        a21 = ab(2) - f2 * ab(3)
        a22 = ac(2) - f2 * ac(3)
      ELSE
        f1 = 0
        f2 = 0 
        b1 = 0
        b2 = 0
        a11 = 0
        a12 = 0
        a21 = 0
        a22 = 0
      END IF
c
      det = a11*a22 - a12*a21
      IF(det .NE. zero) THEN
        r   = (a22 * b1 - a12 * b2) / det
        s   = (a11 * b2 - a21 * b1) / det
      ELSE
        r = 0
        s = 0
      END IF
      ! P = projection of axe origin N1 on ABC plane following direction DIR
      p(1)= a(1) + r*ab(1) + s*ac(1)
      p(2)= a(2) + r*ab(2) + s*ac(2)
      p(3)= a(3) + r*ab(3) + s*ac(3)
      ap(1) = p(1) - a(1)
      ap(2) = p(2) - a(2)
      ap(3) = p(3) - a(3)
      bp(1) = p(1) - b(1)
      bp(2) = p(2) - b(2)
      bp(3) = p(3) - b(3)
      cp(1) = p(1) - c(1)
      cp(2) = p(2) - c(2)
      cp(3) = p(3) - c(3)
      r1 = sqrt(ap(1)**2 + ap(2)**2 + ap(3)**2)
      r2 = sqrt(bp(1)**2 + bp(2)**2 + bp(3)**2)
      r3 = sqrt(cp(1)**2 + cp(2)**2 + cp(3)**2)
      cosp = (dir(1) * ap(1) + dir(2) * ap(2) + dir(3) * ap(3)) / max(r1,em20)
      sinp = max(sqrt(one - cosp**2), em20)
      p1   = zero
      p2   = zero
      p3   = zero
      xx(2) = tt / xfact     ! time
c
c     pressure in 1st point of 3N sub-segment
      cosp = (dir(1) * ap(1) + dir(2) * ap(2) + dir(3) * ap(3)) / max(r1,em20)
      sinp = max(sqrt(one - cosp**2), em20)
      r1   = r1 * sinp
      IF (r1 <=  rmax) THEN
        xx(1) = r1/xfacr
        CALL table_interp_dydx(table(ifunc),xx,2,p1,dydx)
      END IF
c     pressure in 2nd point of 3N sub-segment
      cosp = (dir(1) * bp(1) + dir(2) * bp(2) + dir(3) * bp(3)) / max(r2,em20)
      sinp = max(sqrt(one - cosp**2), em20)
      r2   = r2 * sinp
      IF (r2 <= rmax) THEN
        xx(1) = r2/xfacr
        CALL table_interp_dydx(table(ifunc),xx,2,p2,dydx)
      END IF
c     pressure in 3rd point of 3N sub-segment
      cosp = (dir(1) * cp(1) + dir(2) * cp(2) + dir(3) * cp(3)) / max(r3,em20)
      sinp = max(sqrt(one - cosp**2), em20)
      r3   = r3 * sinp
      IF (r3 <= rmax) THEN
        xx(1) = r3/xfacr
        CALL table_interp_dydx(table(ifunc),xx,2,p3,dydx)
      END IF
      press = (p1 + p2 + p3) * third
      IF (press > zero) THEN
c       calculate cos(angle) between cylinder axis and segment normal to scale pressure
        nx = ab(2) * ac(3) - ab(3) *  ac(2)
        ny = ab(3) * ac(1) - ab(1) *  ac(3)
        nz = ab(1) * ac(2) - ab(2) *  ac(1)
        norm  = sqrt(nx**2 + ny**2 + nz**2)
        cosp  = abs(nx * alpha + ny * beta + nz * gamma) / norm
        press = press * cosp
      END IF
c-----------      
      RETURN
      END