RFLU_ModInCellTest.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 in-cell test.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModInCellTest.F90,v 1.4 2008/12/06 08:44:22 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modincelltest

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

  IMPLICIT NONE
   
  PRIVATE
  PUBLIC :: rflu_ict_computetolerance, & 
            rflu_ict_testfaceplanar, & 
            rflu_ict_testfacequadbilinear, & 
            rflu_ict_testincell, & 
            rflu_ict_testincellfancy, & 
            rflu_ict_testincelllohner
      
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
     
  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: RFLU_ModInCellTest.F90,v $ $Revision: 1.4 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  




! *****************************************************************************
!
! Purpose: Compute tolerance for quadrilateral faces for in-cell test. 
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region data
!
! Output: None.
!
! Notes: None. 
!
! *****************************************************************************

SUBROUTINE rflu_ict_computetolerance(pRegion)
   
  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: ifg,ifl,ipatch,ivg,ivl
  REAL(RFREAL) :: eps,epsmax,nx,ny,nz,safetyfactor,xc,yc,zc
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_ComputeTolerance',&
  'RFLU_ModInCellTest.F90')

  IF ( global%myProcid == masterproc .AND. &
       global%verbLevel > verbose_low ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name, & 
                             'Computing tolerance for in-cell test...'
    WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', & 
                                     pregion%iRegionGlobal
  END IF ! global%myProcid

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

  pgrid => pregion%grid

  epsmax = -huge(1.0_rfreal)

  safetyfactor = 2.0_rfreal

! *****************************************************************************
! Compute plane equation defect 
! *****************************************************************************

  DO ifg = 1,pgrid%nFaces
    IF ( pgrid%f2v(4,ifg) /= vert_none ) THEN 
      xc = pgrid%fc(xcoord,ifg)
      yc = pgrid%fc(ycoord,ifg)
      zc = pgrid%fc(zcoord,ifg)
      
      nx = pgrid%fn(xcoord,ifg)
      ny = pgrid%fn(ycoord,ifg)
      nz = pgrid%fn(zcoord,ifg)

      DO ivl = 1,4
        ivg = pgrid%f2v(ivl,ifg)

        eps = (pgrid%xyz(xcoord,ivg) - xc)*nx & 
            + (pgrid%xyz(ycoord,ivg) - yc)*ny & 
            + (pgrid%xyz(zcoord,ivg) - zc)*nz

        epsmax = max(eps,epsmax)
      END DO ! ivl
    END IF ! pGrid%f2v
  END DO ! ifg

  DO ipatch = 1,pgrid%nPatches
    ppatch => pregion%patches(ipatch)

    DO ifl = 1,ppatch%nBFaces
      IF ( ppatch%bf2v(4,ifl) /= vert_none ) THEN 
        xc = ppatch%fc(xcoord,ifl)
        yc = ppatch%fc(ycoord,ifl)
        zc = ppatch%fc(zcoord,ifl)
      
        nx = ppatch%fn(xcoord,ifl)
        ny = ppatch%fn(ycoord,ifl)
        nz = ppatch%fn(zcoord,ifl)

        DO ivl = 1,4
          ivg = ppatch%bv(ppatch%bf2v(ivl,ifl))
  
          eps = (pgrid%xyz(xcoord,ivg) - xc)*nx & 
              + (pgrid%xyz(ycoord,ivg) - yc)*ny & 
              + (pgrid%xyz(zcoord,ivg) - zc)*nz

          epsmax = max(eps,epsmax)
        END DO ! ivl
      END IF ! pPatch%bf2v
    END DO ! ifl
  END DO ! iPatch

! *****************************************************************************
! Write info
! *****************************************************************************

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel > verbose_low ) THEN
    WRITE(stdout,'(A,3X,A,1X,E13.6)') solver_name, & 
                                      'Maximum face planarity defect:',epsmax 
  END IF ! global%myProcid

! *****************************************************************************
! Set new tolerance
! *****************************************************************************

  IF ( epsmax > pregion%mixtInput%tolerICT ) THEN 
    pregion%mixtInput%tolerICT = safetyfactor*epsmax

    IF ( global%myProcid == masterproc .AND. & 
         global%verbLevel > verbose_low ) THEN
      WRITE(stdout,'(A,3X,A,12X,E13.6)') solver_name,'Tolerance reset to:', &
                                         pregion%mixtInput%tolerICT
    END IF ! global%myProcid
  END IF ! epsMax 

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

  IF ( global%myProcid == masterproc .AND. &
       global%verbLevel > verbose_low ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name, & 
                             'Computing tolerance for in-cell test done.'
  END IF ! global%myProcid

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ict_computetolerance







! *****************************************************************************
!
! Purpose: Compute contribution to in-cell test from planar triangular or
!   quadrilateral face. 
!
! Description: Compute dot product of face normal vector and relative vector
!   between face centroid and location in question.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   icg                 Global cell index
!   iPatch              Patch index
!   ifg                 Face index
!
! Output: 
!   dotp                Dot product
!
! Notes: 
!   1. Only applicable to planar faces.
!
! *****************************************************************************

SUBROUTINE rflu_ict_testfaceplanar(pRegion,xLoc,yLoc,zLoc,icg,iPatch,ifg,dotp)
   
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: icg,ifg,ipatch
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  REAL(RFREAL), INTENT(OUT) :: dotp
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2
  REAL(RFREAL) :: fnx,fny,fnz,xcofg,ycofg,zcofg
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_TestFacePlanar',&
  'RFLU_ModInCellTest.F90')

! =============================================================================
! Set grid pointer
! =============================================================================

  pgrid => pregion%grid

! *****************************************************************************
! Compute dot product
! *****************************************************************************

  IF ( ipatch == 0 ) THEN ! interior face
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

    xcofg = pgrid%fc(xcoord,ifg)
    ycofg = pgrid%fc(ycoord,ifg)
    zcofg = pgrid%fc(zcoord,ifg)                    

    IF ( c1 == icg ) THEN 
      fnx =  pgrid%fn(xcoord,ifg)
      fny =  pgrid%fn(ycoord,ifg)
      fnz =  pgrid%fn(zcoord,ifg)
    ELSE IF ( c2 == icg ) THEN
      fnx = -pgrid%fn(xcoord,ifg)
      fny = -pgrid%fn(ycoord,ifg)
      fnz = -pgrid%fn(zcoord,ifg)
    ELSE ! defensive programming
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! c1                    
  ELSE ! boundary face  
    ppatch => pregion%patches(ipatch)

    xcofg = ppatch%fc(xcoord,ifg)
    ycofg = ppatch%fc(ycoord,ifg)
    zcofg = ppatch%fc(zcoord,ifg)                    

    fnx = ppatch%fn(xcoord,ifg)
    fny = ppatch%fn(ycoord,ifg)
    fnz = ppatch%fn(zcoord,ifg)          
  END IF ! iPatch

  dotp = (xcofg-xloc)*fnx + (ycofg-yloc)*fny + (zcofg-zloc)*fnz

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ict_testfaceplanar






! ******************************************************************************
!
! Purpose: Compute contribution to in-cell test from bilinear quadrilateral 
!   face. 
!
! Description: A point is located inside a cell if the dot products of the 
!   outward-facing face normal vectors and the relative vector of the face-
!   centroids and the point in question are all positive.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   iPatch              Patch index
!   ifg                 Face index
!
! Output: 
!   dotp                Dot product indicating whether point inside cell
!
! Notes: 
!   1. Return from this routine early if detect negative dot product.
!
! ******************************************************************************

SUBROUTINE rflu_ict_testfacequadbilinear(pRegion,xLoc,yLoc,zLoc,iPatch,ifg,dotp)
      
  USE rflu_modbilinearpatch    
      
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: ifg,ipatch
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  REAL(RFREAL), INTENT(OUT) :: dotp
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: v1,v2,v3,v4
  REAL(RFREAL) :: ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,nx,ny,nz,u,v,x,x1,x2, & 
                  x3,x4,y,y1,y2,y3,y4,z,z1,z2,z3,z4
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_TestFaceQuadBilinear',&
  'RFLU_ModInCellTest.F90')

! ******************************************************************************
! Set pointer
! ******************************************************************************

  pgrid => pregion%grid
  
! ******************************************************************************
! Get vertices
! ******************************************************************************

  IF ( ipatch == 0 ) THEN ! interior face
    v1 = pgrid%f2v(1,ifg)
    v2 = pgrid%f2v(2,ifg)
    v3 = pgrid%f2v(3,ifg)
    v4 = pgrid%f2v(4,ifg)    

    IF ( v4 == vert_none ) THEN ! defensive coding
! TEMPORARY    
      WRITE(*,*) 'ERROR!'
      stop
! END TEMPORARY
    END IF ! v4                                   
  ELSE ! boundary face  
    ppatch => pregion%patches(ipatch)

    v1 = ppatch%bv(ppatch%bf2v(1,ifg))
    v2 = ppatch%bv(ppatch%bf2v(2,ifg))
    v3 = ppatch%bv(ppatch%bf2v(3,ifg))
                     
    IF ( ppatch%bf2v(4,ifg) == vert_none ) THEN ! defensive coding
! TEMPORARY    
      WRITE(*,*) 'ERROR!'
      stop
! END TEMPORARY
    END IF ! pPatch%bf2v(4,ifg)
    
    v4 = ppatch%bv(ppatch%bf2v(4,ifg))              
  END IF ! iloc

  x1 = pgrid%xyz(xcoord,v1)
  x2 = pgrid%xyz(xcoord,v2)
  x3 = pgrid%xyz(xcoord,v3)   
  x4 = pgrid%xyz(xcoord,v4) 

  y1 = pgrid%xyz(ycoord,v1)
  y2 = pgrid%xyz(ycoord,v2)
  y3 = pgrid%xyz(ycoord,v3)   
  y4 = pgrid%xyz(ycoord,v4) 

  z1 = pgrid%xyz(zcoord,v1)
  z2 = pgrid%xyz(zcoord,v2)
  z3 = pgrid%xyz(zcoord,v3)   
  z4 = pgrid%xyz(zcoord,v4) 

! ******************************************************************************
! Compute geometric terms of parametric patch representation
! ******************************************************************************

  ax = x1 - x2 - x4 + x3
  ay = y1 - y2 - y4 + y3
  az = z1 - z2 - z4 + z3 

  bx = x2 - x1
  by = y2 - y1
  bz = z2 - z1  

  cx = x4 - x1
  cy = y4 - y1
  cz = z4 - z1  

  dx = x1
  dy = y1
  dz = z1
  
! ******************************************************************************
! Compute closest point on bilinear patch 
! ******************************************************************************

  CALL rflu_blin_findclosestpoint(global,ax,ay,az,bx,by,bz,cx,cy,cz, & 
                                  dx,dy,dz,xloc,yloc,zloc,u,v,x,y,z)  
   
! ******************************************************************************
! Compute intersection distance (here referred to simply as dot product)
! ******************************************************************************

  IF ( (u < 0.0_rfreal .OR. u > 1.0_rfreal) .OR. & 
       (v < 0.0_rfreal .OR. v > 1.0_rfreal) ) THEN 
    dotp = crazy_value_int ! Must be large negative number
  ELSE 
    CALL rflu_blin_computenormal(global,ax,ay,az,bx,by,bz,cx,cy,cz, &
                                 dx,dy,dz,u,v,nx,ny,nz)  

    dotp = (x-xloc)*nx + (y-yloc)*ny + (z-zloc)*nz
  END IF ! u

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ict_testfacequadbilinear






! ******************************************************************************
!
! Purpose: Determine whether point is inside cell.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   icg                 Global cell index
!
! Output: 
!   RFLU_ICT_TestInCell = .TRUE.     if point inside cell
!   RFLU_ICT_TestInCell = .FALSE.    if point outside cell
!
! Notes: 
!   1. Return from this routine early if detect negative dot product.
!
! ******************************************************************************

LOGICAL FUNCTION rflu_ict_testincell(pRegion,xLoc,yLoc,zLoc,icg)
   
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: icg
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: cntr,icl,ict,ifg,ifl,ipatch,nfaces,v4
  REAL(RFREAL) :: dotp,toler
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_TestInCell',&
  'RFLU_ModInCellTest.F90')

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

  pgrid => pregion%grid
  
  rflu_ict_testincell = .false.

  toler = -pregion%mixtInput%tolerICT 

! ******************************************************************************
! Check whether point is inside cell
! ******************************************************************************

  ict = pgrid%cellGlob2Loc(1,icg) ! cell type
  icl = pgrid%cellGlob2Loc(2,icg) ! local cell index

  cntr = 0

  SELECT CASE ( ict )

! ==============================================================================
!   Tetrahedron
! ==============================================================================

    CASE ( cell_type_tet )
      nfaces = 4
    
      DO ifl = 1,nfaces
        ipatch = pgrid%tet2f(1,ifl,icl)
        ifg    = pgrid%tet2f(2,ifl,icl)        

        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)

        IF ( dotp >= toler ) THEN 
          cntr = cntr + 1
        ELSE 
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifc     

! ==============================================================================
!   Hexahedron
! ==============================================================================

    CASE ( cell_type_hex ) 
      nfaces = 6
    
      DO ifl = 1,nfaces
        ipatch = pgrid%hex2f(1,ifl,icl)
        ifg    = pgrid%hex2f(2,ifl,icl)        

! TEMPORARY
        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)
!        CALL RFLU_ICT_TestFaceQuadBilinear(pRegion,xLoc,yLoc,zLoc,iPatch,ifg, &
!                                           dotp)
! END TEMPORARY

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifc   

! ==============================================================================
!   Prism
! ==============================================================================

    CASE ( cell_type_pri )                 
      nfaces = 5
    
      DO ifl = 1,nfaces
        ipatch = pgrid%pri2f(1,ifl,icl)
        ifg    = pgrid%pri2f(2,ifl,icl)        

        IF ( ipatch == 0 ) THEN 
          v4 = pgrid%f2v(4,ifg)
        ELSE IF ( ipatch > 0 ) THEN
          ppatch => pregion%patches(ipatch)
         
          v4 = ppatch%bf2v(4,ifg)
        ELSE 
          CALL errorstop(global,err_reached_default,__line__)
        END IF ! iPatch        

        IF ( v4 == vert_none ) THEN 
          CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg, &
                                       dotp)  
        ELSE
! TEMPORARY
          CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg, &
                                       dotp)  
!          CALL RFLU_ICT_TestFaceQuadBilinear(pRegion,xLoc,yLoc,zLoc,iPatch, &
!                                             ifg,dotp)
! END TEMPORARY
        END IF ! v4

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifc 

! ==============================================================================
!   Pyramid
! ==============================================================================

    CASE ( cell_type_pyr ) 
      nfaces = 5
    
      DO ifl = 1,nfaces
        ipatch = pgrid%pyr2f(1,ifl,icl)
        ifg    = pgrid%pyr2f(2,ifl,icl)        

        IF ( ipatch == 0 ) THEN 
          v4 = pgrid%f2v(4,ifg)
        ELSE IF ( ipatch > 0 ) THEN
          ppatch => pregion%patches(ipatch)
         
          v4 = ppatch%bf2v(4,ifg)
        ELSE 
          CALL errorstop(global,err_reached_default,__line__)
        END IF ! iPatch        

        IF ( v4 == vert_none ) THEN 
          CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg, &
                                       dotp)  
        ELSE
! TEMPORARY
          CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg, &
                                       dotp)  
!          CALL RFLU_ICT_TestFaceQuadBilinear(pRegion,xLoc,yLoc,zLoc,iPatch, &
!                                             ifg,dotp)
! END TEMPORARY
        END IF ! v4
        
        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifc 

! ==============================================================================
!   Default
! ==============================================================================

    CASE default
      CALL errorstop(global,err_reached_default,__line__)                      
  END SELECT ! ict        

! ==============================================================================
! Set RFLU_ICT_TestInCell to TRUE if have nFaces positive dot products
! ==============================================================================

  IF ( cntr == nfaces ) THEN 
    rflu_ict_testincell = .true. 
  END IF ! cntr

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

  CALL deregisterfunction(global)

END FUNCTION rflu_ict_testincell







! *****************************************************************************
!
! Purpose: Determine whether point is inside cell and return information about
!   which face which might have lead to failure.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   icg                 Global cell index
!
! Output: 
!   testInCell          True if location in cell, false otherwise
!   iPatchOut           Patch index of face which failed test
!   ifgOut              Global index of face which failed test
!
! Notes: 
!   1. Return from this routine early if detect negative dot product.
!
! *****************************************************************************

SUBROUTINE rflu_ict_testincellfancy(pRegion,xLoc,yLoc,zLoc,icg,testInCell, & 
                                    ipatchout,ifgout)
      
  IMPLICIT NONE

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

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

  LOGICAL, INTENT(OUT) :: testincell
  INTEGER, INTENT(IN) :: icg
  INTEGER, INTENT(OUT) :: ipatchout,ifgout
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: cntr,icl,ict,ifg,ifl,ipatch,nfaces
  REAL(RFREAL) :: dotp,toler
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_TestInCellFancy',&
  'RFLU_ModInCellTest.F90')

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

  pgrid => pregion%grid
  
  testincell = .false.

  toler = -pregion%mixtInput%tolerICT

! *****************************************************************************
! Check whether point is inside cell
! *****************************************************************************

  ict = pgrid%cellGlob2Loc(1,icg) ! cell type
  icl = pgrid%cellGlob2Loc(2,icg) ! local cell index

  cntr = 0

  SELECT CASE ( ict )
    CASE ( cell_type_tet )
      nfaces = 4
    
      DO ifl = 1,nfaces
        ipatch = pgrid%tet2f(1,ifl,icl)
        ifg    = pgrid%tet2f(2,ifl,icl)        

        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          ipatchout = ipatch
          ifgout    = ifg
        
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifl    
    CASE ( cell_type_hex ) 
      nfaces = 6
    
      DO ifl = 1,nfaces
        ipatch = pgrid%hex2f(1,ifl,icl)
        ifg    = pgrid%hex2f(2,ifl,icl)        

        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          ipatchout = ipatch
          ifgout    = ifg

          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifl   
    CASE ( cell_type_pri )                 
      nfaces = 5
    
      DO ifl = 1,nfaces
        ipatch = pgrid%pri2f(1,ifl,icl)
        ifg    = pgrid%pri2f(2,ifl,icl)        

        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          ipatchout = ipatch
          ifgout    = ifg
                
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifl 
    CASE ( cell_type_pyr ) 
      nfaces = 5
    
      DO ifl = 1,nfaces
        ipatch = pgrid%pyr2f(1,ifl,icl)
        ifg    = pgrid%pyr2f(2,ifl,icl)        

        CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg,dotp)

        IF ( dotp > toler ) THEN 
          cntr = cntr + 1
        ELSE 
          ipatchout = ipatch
          ifgout    = ifg
                
          CALL deregisterfunction(global)
          RETURN
        END IF ! dotp
      END DO ! ifl 
    CASE default
      CALL errorstop(global,err_reached_default,__line__)                      
  END SELECT          

! =============================================================================
! Set testInCell to TRUE if have nFaces positive dot products
! =============================================================================

  IF ( cntr == nfaces ) THEN 
    testincell = .true. 
    
    ipatchout = crazy_value_int
    ifgout    = crazy_value_int
  END IF ! cntr

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ict_testincellfancy







! ******************************************************************************
!
! Purpose: Determine whether point is inside cell and return information about
!   which face failed the most, so to speak.
!
! Description: None.
!
! Input:
!   pRegion             Pointer to region data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   icg                 Global cell index
!
! Output: 
!   testInCell          True if location in cell, false otherwise
!   ilocOut             Location of face which failed test
!   ifgOut              Global index of face which failed test
!
! Notes: 
!   1. DO NOT return from this routine early if detect negative dot product.
!   2. See R. Lohner and J. Ambrosiano, A Vectorized Particle Tracer for 
!      Unstructured Grids, J. Comp. Phys., Vol. 91, pp. 22-31, 1990.
!
! ******************************************************************************

SUBROUTINE rflu_ict_testincelllohner(pRegion,xLoc,yLoc,zLoc,icg,testInCell, & 
                                     ipatchout,ifgout)
   
  IMPLICIT NONE

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

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

  LOGICAL, INTENT(OUT) :: testincell
  INTEGER, INTENT(IN) :: icg
  INTEGER, INTENT(OUT) :: ipatchout,ifgout
  REAL(RFREAL), INTENT(IN) :: xloc,yloc,zloc
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: cntr,icl,ict,ifg,ifl,ipatch,nfaces
  INTEGER, DIMENSION(1) :: idotpmin
  INTEGER, DIMENSION(:,:), POINTER :: pc2f
  REAL(RFREAL) :: toler
  REAL(RFREAL), DIMENSION(6) :: dotp
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ICT_TestInCellLohner',&
  'RFLU_ModInCellTest.F90')

! ==============================================================================
! Set pointer and initialize 
! ==============================================================================

  pgrid => pregion%grid
  
  testincell = .false.

  toler = -pregion%mixtInput%tolerICT

! ******************************************************************************
! Select cell type and set pointer to cell-to-face connectivity array
! ******************************************************************************

  ict = pgrid%cellGlob2Loc(1,icg) ! cell type
  icl = pgrid%cellGlob2Loc(2,icg) ! local cell index

  SELECT CASE ( ict ) 
    CASE ( cell_type_tet )
      pc2f => pgrid%tet2f(:,:,icl)
    CASE ( cell_type_hex ) 
      pc2f => pgrid%hex2f(:,:,icl)
    CASE ( cell_type_pri ) 
      pc2f => pgrid%pri2f(:,:,icl)      
    CASE ( cell_type_pyr ) 
      pc2f => pgrid%pyr2f(:,:,icl)          
    CASE default 
      CALL errorstop(global,err_reached_default,__line__)      
  END SELECT ! ict
  
  nfaces = SIZE(pc2f,2)

! ******************************************************************************
! Check whether point is inside cell
! ******************************************************************************

  cntr = 0

! ==============================================================================
! Loop over faces of cell. NOTE need to initialize dotp for non-existent faces.
! ==============================================================================
    
  DO ifl = 1,nfaces
    ipatch = pc2f(1,ifl)
    ifg  = pc2f(2,ifl)        

    CALL rflu_ict_testfaceplanar(pregion,xloc,yloc,zloc,icg,ipatch,ifg, &
                                 dotp(ifl))

    IF ( dotp(ifl) > toler ) THEN 
      cntr = cntr + 1
    END IF ! dotp
  END DO ! ifl   

  IF ( nfaces < 6 ) THEN 
    DO ifl = nfaces+1,6
      dotp(ifl) = huge(1.0_rfreal) 
    END DO ! ifl
  END IF ! nFaces

! ==============================================================================
! Set testInCell to TRUE if have nFaces positive dot products
! ==============================================================================

  IF ( cntr == nfaces ) THEN 
    testincell = .true. 
    
    ipatchout = crazy_value_int
    ifgout    = crazy_value_int
  ELSE 
    idotpmin = minloc(dotp(:))
  
    ipatchout = pc2f(1,idotpmin(1))
    ifgout    = pc2f(2,idotpmin(1))  
  END IF ! cntr

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_ict_testincelllohner





! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modincelltest


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModInCellTest.F90,v $
! Revision 1.4  2008/12/06 08:44:22  mtcampbe
! Updated license.
!
! Revision 1.3  2008/11/19 22:17:33  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.2  2006/04/07 15:19:19  haselbac
! Removed tabs
!
! Revision 1.1  2005/12/24 21:17:50  haselbac
! Initial revision
!
! ******************************************************************************