press_seg3.F File Reference

#include "implicit_f.inc"
#include "com04_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	press_seg3 (a, b, c, n1, dir, ifunc, table, xfacr, xfact, press)

Function/Subroutine Documentation

press_seg3()

subroutine press_seg3	(	intent(in)	a,
		intent(in)	b,
		intent(in)	c,
		intent(in)	n1,
		intent(in)	dir,
		integer, intent(in)	ifunc,
		type (ttable), dimension(ntable), intent(in)	table,
		intent(in)	xfacr,
		intent(in)	xfact,
		intent(out)	press )

Definition at line 35 of file press_seg3.F.

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