RFLU_ModBilinearPatch.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * 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.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Suite of routines related to bilinear patches.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModBilinearPatch.F90,v 1.5 2008/12/06 08:44:20 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modbilinearpatch

  USE modglobal, ONLY: t_global 
  USE modparameters
  USE moddatatypes
  USE moderror
  USE modbndpatch, ONLY: t_patch
  USE moddatastruct, ONLY: t_region
  USE modgrid, ONLY: t_grid
  USE modmpi

  IMPLICIT NONE
   
  PRIVATE
  PUBLIC :: rflu_blin_computenormal, &
            rflu_blin_computexsectline, &  
            rflu_blin_findclosestpoint
      
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
     
  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: RFLU_ModBilinearPatch.F90,v $ $Revision: 1.5 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  





! ******************************************************************************
!
! Purpose: Compute normal vector of bilinear patch at given (u,v) coordinates.
!
! Description: None.
!
! Input:
!   global      Pointer to global data
!   ax,ay,az    Components of bilinear patch parametric representation
!   bx,by,bz    Components of bilinear patch parametric representation
!   cx,cy,cz    Components of bilinear patch parametric representation
!   dx,dy,dz    Components of bilinear patch parametric representation
!   u           u-coordinate
!   v           v-coordinate
!
! Output: 
!   nx          x-component of normal vector
!   ny          y-component of normal vector
!   nz          z-component of normal vector
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_blin_computenormal(global,ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz, &
                                   u,v,nx,ny,nz)

  IMPLICIT NONE

! ******************************************************************************
! Declarations and definitions
! ******************************************************************************

! ==============================================================================
! Parameters      
! ==============================================================================

  REAL(RFREAL), INTENT(IN) :: ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,u,v
  REAL(RFREAL), INTENT(OUT) :: nx,ny,nz
  TYPE(t_global), POINTER :: global   

! ==============================================================================
! Locals      
! ==============================================================================

  REAL(RFREAL) :: imag

! ******************************************************************************
! Start
! ******************************************************************************

  CALL registerfunction(global,'RFLU_BLIN_ComputeNormal',&
  'RFLU_ModBilinearPatch.F90')

! ******************************************************************************
! Compute normal
! ******************************************************************************

  nx = v*(ay*cz - az*cy) - u*(ay*bz - az*by) + by*cz - bz*cy
  ny = v*(az*cx - ax*cz) - u*(az*bx - ax*bz) + bz*cx - bx*cz
  nz = v*(ax*cy - ay*cx) - u*(ax*by - ay*bx) + bx*cy - by*cx        

  imag = 1.0_rfreal/sqrt(nx*nx + ny*ny + nz*nz)

  nx = imag*nx
  ny = imag*ny
  nz = imag*nz

! ******************************************************************************
! End
! ******************************************************************************  

  CALL deregisterfunction(global)  

END SUBROUTINE rflu_blin_computenormal










! ******************************************************************************
!
! Purpose: Compute intersection distance of given location along given line 
!   vector and bilinear face. 
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   ex,ey,ez            x-, y-, and z-components of unit line vector
!   iPatch              Face location 
!   ifg                 Face index
!
! Output:
!   nt                  Number of intersections
!   t                   Array of intersection distances
!
! Notes: 
!   1. The line vector MUST be a unit line vector. If that is not correct, the
!      distance between the given point and the intersection of the line with 
!      the faces of the given cell will not be computed correctly.
!   2. For details on derivation, see: S.D. Ramsey, K. Potter, and C. Hansen, 
!      Ray Bilinear Patch Intersection, ACM J. Graphics Tools, Vol. 9, No. 3, 
!      pp. 41-47, 2004. 
!
! ******************************************************************************

SUBROUTINE rflu_blin_computexsectline(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
                                      ipatch,ifg,nt,t)
     
  USE modtools, ONLY: swaprfreals   
     
  IMPLICIT NONE

! ******************************************************************************
! Declarations and definitions
! ******************************************************************************

! ==============================================================================
! Arguments
! ==============================================================================

  INTEGER, INTENT(IN) :: icg,ifg,ipatch
  INTEGER, INTENT(INOUT) :: nt
  REAL(RFREAL), INTENT(IN) :: ex,ey,ez
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc    
  REAL(RFREAL), INTENT(OUT) :: t(2)
  TYPE(t_region), POINTER :: pregion

! ==============================================================================
! Locals
! ==============================================================================

  INTEGER :: iv,nv,v1g,v1l,v2g,v2l,v3g,v3l,v4g,v4l
  REAL(RFREAL) :: a,ae,ax,ay,az,a1,a2,b,be,bx,by,bz,b1,b2,c,ce,cx,cy,cz,c1, &
                  c2,d,de,denom1,denom2,disc,dotp,dx,dy,dz,d1,d2, &
                  fuzzytolerance,px,py,pz,q,re,u,xcofg,x1,x2,x3,x4, &
                  ycofg,y1,y2,y3,y4,zcofg,z1,z2,z3,z4
  REAL(RFREAL), DIMENSION(2) :: v
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

! DEBUG
  INTEGER :: ivg
  REAL(RFREAL) :: nx,ny,nz,a3,b3,c3,d3
! END DEBUG

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global

  CALL registerfunction(global,'RFLU_BLIN_ComputeXSectLine',&
  'RFLU_ModBilinearPatch.F90')

! ==============================================================================
! Set grid pointer and initialize variables
! ==============================================================================

  pgrid => pregion%grid

  nt = 0
  
  fuzzytolerance = -pregion%mixtInput%tolerICT

  t(1:2) = REAL(crazy_value_int,kind=rfreal)

! ******************************************************************************
! Get face vertices and coordinates
! ******************************************************************************

! DEBUG
  WRITE(0,*) '###100',pregion%iRegionGlobal,ipatch,ifg
  WRITE(0,*) '###110',pregion%iRegionGlobal,xloc,yloc,zloc
  WRITE(0,*) '###120',pregion%iRegionGlobal,ex,ey,ez
  
!  DO ivg = 1,pGrid%nVert
!    WRITE(*,'(1X,I2,3(1X,E13.6))') ivg,pGrid%xyz(:,ivg)
!  END DO ! ivg
! END DEBUG  

  IF ( ipatch == 0 ) THEN 
    v1g = pgrid%f2v(1,ifg)
    v2g = pgrid%f2v(2,ifg)
    v3g = pgrid%f2v(3,ifg)
    v4g = pgrid%f2v(4,ifg)  
    
    IF ( v4g == vert_none ) THEN 
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! v4g        
  ELSE IF ( ipatch > 0 ) THEN 
    ppatch => pregion%patches(ipatch)
    
    v1l = ppatch%bf2v(1,ifg)
    v2l = ppatch%bf2v(2,ifg)
    v3l = ppatch%bf2v(3,ifg)
    v4l = ppatch%bf2v(4,ifg) 

    IF ( v4l == vert_none ) THEN 
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! v4l
      
    v1g = ppatch%bv(v1l)
    v2g = ppatch%bv(v2l)
    v3g = ppatch%bv(v3l)
    v4g = ppatch%bv(v4l)                    
  ELSE
    CALL errorstop(global,err_reached_default,__line__)
  END IF ! iPatch
  
! DEBUG
  WRITE(0,*) '###130',pregion%iRegionGlobal,v1g,v2g,v3g,v4g
! END DEBUG  
  
  x1 = pgrid%xyz(xcoord,v1g)
  x2 = pgrid%xyz(xcoord,v2g)
  x3 = pgrid%xyz(xcoord,v3g)
  x4 = pgrid%xyz(xcoord,v4g)      

  y1 = pgrid%xyz(ycoord,v1g)
  y2 = pgrid%xyz(ycoord,v2g)
  y3 = pgrid%xyz(ycoord,v3g)
  y4 = pgrid%xyz(ycoord,v4g) 
  
  z1 = pgrid%xyz(zcoord,v1g)
  z2 = pgrid%xyz(zcoord,v2g)
  z3 = pgrid%xyz(zcoord,v3g)
  z4 = pgrid%xyz(zcoord,v4g)   

! DEBUG
  WRITE(0,*) '@@@131a',pregion%iRegionGlobal,x1,x2,x3,x4
  WRITE(0,*) '@@@131b',pregion%iRegionGlobal,y1,y2,y3,y4
  WRITE(0,*) '@@@131c',pregion%iRegionGlobal,z1,z2,z3,z4
! END DEBUG

! ******************************************************************************
! Compute geometric terms
! ******************************************************************************

  re = xloc*ex + yloc*ey + zloc*ez

  ax = x3 - x2 - x4 + x1
  ay = y3 - y2 - y4 + y1
  az = z3 - z2 - z4 + z1  
  ae = ax*ex + ay*ey + az*ez

  bx = x2 - x1
  by = y2 - y1
  bz = z2 - z1   
  be = bx*ex + by*ey + bz*ez 
  
  cx = x4 - x1
  cy = y4 - y1
  cz = z4 - z1    
  ce = cx*ex + cy*ey + cz*ez

  dx = x1
  dy = y1
  dz = z1        
  de = dx*ex + dy*ey + dz*ez

  a1 = az - ay - ae*(ez - ey)
  a2 = az - ax - ae*(ez - ex)
  a3 = a2 - a1
  
  b1 = bz - by - be*(ez - ey)
  b2 = bz - bx - be*(ez - ex)
  b3 = b2 - b1
  
  c1 = cz - cy - ce*(ez - ey)
  c2 = cz - cx - ce*(ez - ex)
  c3 = c2 - c1
  
  d1 = dz - dy - (zloc - yloc) - (de - re)*(ez - ey)
  d2 = dz - dx - (zloc - xloc) - (de - re)*(ez - ex)
  d3 = d2 - d1

  a = a1*c2 - a2*c1
  b = a1*d2 - a2*d1 + b1*c2 - b2*c1
  c = b1*d2 - b2*d1
  
! DEBUG
  WRITE(0,*) 'a:',ax,ay,az,ae
  WRITE(0,*) 'b:',bx,by,bz,be
  WRITE(0,*) 'c:',cx,cy,cz,ce
  WRITE(0,*) 'd:',dx,dy,dz,de
  WRITE(0,*) 'A-D1:',a1,b1,c1,d1
  WRITE(0,*) 'A-D2:',a2,b2,c2,d2
  WRITE(0,*) 'A-D3:',a3,b3,c3,d3
  WRITE(0,*) 'a-c:',a,b,c
! END DEBUG  
  
! ******************************************************************************
! Compute v 
! ******************************************************************************
  
  IF ( abs(a) > 0.0_rfreal ) THEN ! Quadratic equation   
    disc = b*b - 4.0_rfreal*a*c

! DEBUG
  WRITE(0,*) 'disc=',disc
! END DEBUG  

    IF ( disc > 0.0_rfreal ) THEN
      nv = 2
    
      IF ( abs(c) > 0.0_rfreal ) THEN 
        q = -0.5_rfreal*(b + sign(1.0_rfreal,b)*sqrt(disc))

        v(1) = q/a
        v(2) = c/q
      ELSE 
        v(1) = 0.0_rfreal
        v(2) = -b/a 
      END IF ! nv
    ELSE IF ( disc < 0.0_rfreal ) THEN 
      nv = 0

      v(1) = REAL(CRAZY_VALUE_INT,KIND=RFREAL)  
      v(2) = REAL(CRAZY_VALUE_INT,KIND=RFREAL)    
    ELSE 
      nv = 1

      v(1) = -0.5_rfreal*b/a
      v(2) = REAL(crazy_value_int,kind=rfreal)
    END IF ! disc
  ELSE IF ( abs(b) > 0.0_rfreal ) THEN ! Linear equation
    nv = 1
    
    v(1) = -c/b
  ELSE 
    nv = 0
  END IF ! ABS(a)

! DEBUG
  WRITE(0,*) '###132',pregion%iRegionGlobal,v(:)
! END DEBUG  
  

! ******************************************************************************
! Compute intersection distance
! ******************************************************************************

  DO iv = 1,nv

! ==============================================================================
!   Check value of v, proceed only if within range [0,1]
! ==============================================================================

    IF ( (v(iv) < 0.0_rfreal) .OR. (v(iv) > 1.0_rfreal) ) THEN 
      cycle
    END IF ! v(iv)
    
! ==============================================================================
!   Compute u, proceed only if within range [0,1]
! ==============================================================================

    denom1 = v(iv)*a2 + b2
    denom2 = denom1 - (v(iv)*a1 + b1)
    
    IF ( abs(denom2) > abs(denom1) ) THEN 
      u = (v(iv)*(c1 - c2) + d1 - d2)/denom2
    ELSE 
      u = -(v(iv)*c2 + d2)/denom1
    END IF ! ABS(denom1)
    
! DEBUG
  WRITE(0,*) '###134',pregion%iRegionGlobal,iv,v(iv),u
! END DEBUG  
    
    IF ( (u < 0.0_rfreal) .OR. (u > 1.0_rfreal) ) THEN 
      cycle
    END IF ! u
    
! ==============================================================================
!   Compute intersection point 
! ==============================================================================
   
    px = u*v(iv)*ax + u*bx + v(iv)*cx + dx
    py = u*v(iv)*ay + u*by + v(iv)*cy + dy    
    pz = u*v(iv)*az + u*bz + v(iv)*cz + dz

! DEBUG
  WRITE(0,*) '###135',pregion%iRegionGlobal,px,py,pz
! END DEBUG  
    
! ==============================================================================
!   Compute normal at intersection point and ensure it is an outward normal
! ==============================================================================
    
    CALL rflu_blin_computenormal(global,ax,ay,az,bx,by,bz,cx,cy,cz, &
                                 dx,dy,dz,u,v(iv),nx,ny,nz)
  
    xcofg = pgrid%cofg(xcoord,icg)
    ycofg = pgrid%cofg(ycoord,icg)
    zcofg = pgrid%cofg(zcoord,icg)  
  
    dotp = (px - xcofg)*nx + (py - ycofg)*ny + (pz - zcofg)*nz
  
! DEBUG
  WRITE(0,*) '###136',pregion%iRegionGlobal,nx,ny,nz,dotp
! END DEBUG  
    
    IF ( dotp < 0.0_rfreal ) THEN
      nx = -nx
      ny = -ny
      nz = -nz
    END IF ! dotp
    
! ==============================================================================
!   Compute intersection distance. NOTE expression for intersection
!   distance differs from that in the paper - it avoids IFs and divisions.
! ==============================================================================
  
    dotp = (px - xloc)*nx + (py - yloc)*ny + (pz - zloc)*nz
  
! DEBUG
  WRITE(0,*) '###137',pregion%iRegionGlobal,dotp
! END DEBUG  
    
    IF ( dotp > fuzzytolerance ) THEN ! Inside cell 
! DEBUG
  WRITE(0,*) '###138',pregion%iRegionGlobal,ex*nx + ey*ny + ez*nz
! END DEBUG  
    
      IF ( ex*nx + ey*ny + ez*nz > 0.0_rfreal ) THEN ! In direction of traj
        nt = nt + 1
    
        t(nt) = (px - xloc)*ex + (py - yloc)*ey + (pz - zloc)*ez
      END IF ! ex*nx
    END IF ! dotp
  END DO ! iv
  
! ******************************************************************************
! Swap so first entry is that with smaller intersection distance
! ******************************************************************************

  IF ( nt > 0 ) THEN 
    IF ( abs(t(2)) < abs(t(1)) ) THEN 
      CALL swaprfreals(t(1),t(2))
    END IF ! ABS(t(2)) 
  END IF ! nt

! DEBUG
  WRITE(0,*) '###140',pregion%iRegionGlobal,nt,t(1:2)
! END DEBUG  
  
! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_blin_computexsectline








! ******************************************************************************
!
! Purpose: Given a point, compute closest point on bilinear patch.
!
! Description: None.
!
! Input:
!   global      Pointer to global data
!   ax,ay,az    Components of bilinear patch parametric representation
!   bx,by,bz    Components of bilinear patch parametric representation
!   cx,cy,cz    Components of bilinear patch parametric representation
!   dx,dy,dz    Components of bilinear patch parametric representation
!   xq,yq,zq    Coordinates of point
!
! Output: 
!   u,v         Parametric coordinates of closest point on bilinear patch
!   x,y,z       Coordinates of closest point on bilinear patch
!
! Notes:
!   1. Use Newton-Raphson method to solve non-linear system for (u,v). 
!   2. This routine does NOT check that the closest point lies in the patch, 
!      i.e., that the (u,v) coordinates are in the range [0,1].
!
! ******************************************************************************

  SUBROUTINE rflu_blin_findclosestpoint(global,ax,ay,az,bx,by,bz,cx,cy,cz, & 
                                        dx,dy,dz,xq,yq,zq,u,v,x,y,z)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Parameters      
! ==============================================================================

    REAL(RFREAL), INTENT(OUT) :: u,v,x,y,z
    REAL(RFREAL), INTENT(IN) :: ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,xq,yq,zq
    TYPE(t_global), POINTER :: global   
    
! ==============================================================================
!   Locals      
! ==============================================================================

    INTEGER :: loopcounter
    REAL(RFREAL) :: dfdu,dfdv,dgdu,dgdv,du,dv,dxdu,dxdv,dydu,dydv,dzdu,dzdv, &
                    d2,d2xduv,d2yduv,d2zduv,f,g,idet,tol
    
! ******************************************************************************
!   Start
! ******************************************************************************

    CALL registerfunction(global,'RFLU_BLIN_FindClosestPoint',&
  'RFLU_ModBilinearPatch.F90')


! ******************************************************************************
!   Initialize
! ******************************************************************************

    loopcounter = 0
    
    tol = 1.0e-6_rfreal

    u = 0.5_rfreal ! Initial guess
    v = 0.5_rfreal

! ******************************************************************************
!   Newton-Raphson iteration to compute closest point
! ******************************************************************************
        
    emptyloop: DO 
      loopcounter = loopcounter + 1
    
      x = u*v*ax + u*bx + v*cx + dx 
      y = u*v*ay + u*by + v*cy + dy
      z = u*v*az + u*bz + v*cz + dz            
    
      dxdu = v*ax + bx 
      dydu = v*ay + by
      dzdu = v*az + bz  
      
      dxdv = u*ax + cx
      dydv = u*ay + cy
      dzdv = u*az + cz                  
    
      d2xduv = ax
      d2yduv = ay
      d2zduv = az
    
      d2 = (x-xq)*(x-xq) + (y-yq)*(y-yq) + (z-zq)*(z-zq) 
    
      f = (x-xq)*dxdu + (y-yq)*dydu + (z-zq)*dzdu
      g = (x-xq)*dxdv + (y-yq)*dydv + (z-zq)*dzdv
      
      dfdu = dxdu*dxdu + dydu*dydu + dzdu*dzdu
      dfdv = dxdu*dxdv + dydu*dydv + dzdu*dzdv & 
           + (x-xq)*d2xduv + (y-yq)*d2yduv + (z-zq)*d2zduv 

      dgdu = dfdv 
      dgdv = dxdv*dxdv + dydv*dydv + dzdv*dzdv
      
      idet = 1.0_rfreal/(dfdu*dgdv - dfdv*dgdu)
      
      du = -idet*(dgdv*f - dgdu*g)
      dv = -idet*(dfdu*g - dfdv*f)

      IF ( abs(du) < tol .AND. abs(dv) < tol ) THEN 
        EXIT emptyloop
      ELSE IF ( loopcounter > limit_infinite_loop ) THEN 
        CALL errorstop(global,err_infinite_loop,__line__)
      END IF ! ABS(du) 

      u = u + du 
      v = v + dv      
    END DO emptyloop 

! ******************************************************************************
!   End
! ******************************************************************************  

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_blin_findclosestpoint






! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modbilinearpatch


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModBilinearPatch.F90,v $
! Revision 1.5  2008/12/06 08:44:20  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:31  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2006/04/07 15:19:18  haselbac
! Removed tabs
!
! Revision 1.2  2006/03/25 21:49:14  haselbac
! Renamed SwapFloats to SwapRFREALs
!
! Revision 1.1  2005/12/24 21:17:50  haselbac
! Initial revision
!
! ******************************************************************************