RFLU_ModGeometryTools.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 geometry tools.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModGeometryTools.F90,v 1.5 2008/12/06 08:44:22 mtcampbe Exp $
!
! Copyright: (c) 2005-2007 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modgeometrytools

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

  IMPLICIT NONE
   
  PRIVATE
  PUBLIC :: rflu_computelinecellxsectfast, & 
            rflu_computelinecellxsectsafe, & 
            rflu_testinboundbox, &
            rflu_testvectorcartaxisaligned
      
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
     
  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: RFLU_ModGeometryTools.F90,v $ $Revision: 1.5 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  




! ******************************************************************************
!
! Purpose: Compute intersection of given line vector and faces of given cell
!   and distance between given location and intersection using fast algorithm. 
!
! 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
!   icg                 Global cell index
!
! Output:
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of intersection
!   distOut             Distance from location to intersection 
!   iPatchOut           Face location 
!   ifgOut              Face index 
!
! 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. Distance might be zero if vertices lie on path, so need to include this
!      case in IF statement on dist after calls to routine which computes the
!      distance.
!   3. This is called a fast algorithm because it CANNOT detect and correct 
!      inconsistent input data, which makes it faster than the safe version.
!
! ******************************************************************************

SUBROUTINE rflu_computelinecellxsectfast(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
                                         distout,ipatchout,ifgout)
     
  USE rflu_modbilinearpatch, ONLY: rflu_blin_computexsectline  
     
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: icg
  INTEGER, INTENT(OUT) :: ifgout,ipatchout
  REAL(RFREAL), INTENT(IN) :: ex,ey,ez
  REAL(RFREAL), INTENT(INOUT) :: xloc,yloc,zloc    
  REAL(RFREAL), INTENT(OUT) :: distout
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,icl,ict,ifl,iflout,ifg,ipatch,nfaces,nxsect
  INTEGER, DIMENSION(:,:), POINTER :: pc2f
  REAL(RFREAL) :: denom,dist,fnx,fny,fnz,numer,toler,xcofg,ycofg,zcofg
  REAL(RFREAL), DIMENSION(2) :: xsd  
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ComputeLineCellXSectFast',&
  'RFLU_ModGeometryTools.F90')

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

  pgrid => pregion%grid

  toler = -pregion%mixtInput%tolerICT ! Must be consistent with ICT tolerance

  distout = huge(1.0_rfreal)
  iflout  = crazy_value_int

! ******************************************************************************
! 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)

! ******************************************************************************
! Loop over faces of cell
! ******************************************************************************
  
  DO ifl = 1,nfaces
    ipatch = pc2f(1,ifl)
    ifg    = pc2f(2,ifl)

! ==============================================================================
!   Interior face
! ==============================================================================

    IF ( ipatch == 0 ) THEN 

! ------------------------------------------------------------------------------
!     Triangular face
! ------------------------------------------------------------------------------    
    
! TEMPORARY
!      IF ( pGrid%f2v(4,ifg) == VERT_NONE ) THEN 
! END TEMPORARY
        c1 = pgrid%f2c(1,ifg)
        c2 = pgrid%f2c(2,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                 

        denom = ex*fnx + ey*fny + ez*fnz

        IF ( denom > 0.0_rfreal ) THEN 
          xcofg = pgrid%fc(xcoord,ifg)
          ycofg = pgrid%fc(ycoord,ifg)
          zcofg = pgrid%fc(zcoord,ifg) 

          numer = (xcofg-xloc)*fnx + (ycofg-yloc)*fny + (zcofg-zloc)*fnz      
          numer = max(numer,0.0_rfreal) 

          dist = numer/denom

          IF ( dist < distout ) THEN 
            distout = dist
            iflout  = ifl
          END IF ! dist          
        END IF ! denom 

! ------------------------------------------------------------------------------
!     Quadrilateral face
! ------------------------------------------------------------------------------    
    
! TEMPORARY
!      ELSE 
!        CALL RFLU_BLIN_ComputeXSectLine(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
!                                        iPatch,ifg,nXSect,xsd)        
!                              
!        IF ( nXSect > 0 ) THEN 
!          dist = xsd(1)
!          
!          IF ( dist < distOut ) THEN 
!            distOut = dist
!            iflOut  = ifl
!          END IF ! dist          
!        END IF ! nXSect                                   
!      END IF ! pGrid%f2v     
! END TEMPORARY

! ==============================================================================
!   Boundary face
! ==============================================================================
         
    ELSE IF ( ipatch > 0 ) THEN  
      ppatch => pregion%patches(ipatch)

! ------------------------------------------------------------------------------
!     Triangular face
! ------------------------------------------------------------------------------    

! TEMPORARY
!      IF ( pPatch%bf2v(4,ifg) == VERT_NONE ) THEN 
! END TEMPORARY
        fnx = ppatch%fn(xcoord,ifg)
        fny = ppatch%fn(ycoord,ifg)
        fnz = ppatch%fn(zcoord,ifg)          

        denom = ex*fnx + ey*fny + ez*fnz

        IF ( denom > 0.0_rfreal ) THEN 
          xcofg = ppatch%fc(xcoord,ifg)
          ycofg = ppatch%fc(ycoord,ifg)
          zcofg = ppatch%fc(zcoord,ifg)

          numer = (xcofg-xloc)*fnx + (ycofg-yloc)*fny + (zcofg-zloc)*fnz      
          numer = max(numer,0.0_rfreal) 

          dist = numer/denom

          IF ( dist < distout ) THEN 
            distout = dist
            iflout  = ifl
          END IF ! dist 
        END IF ! denom 

! ------------------------------------------------------------------------------
!     Quadrilateral face
! ------------------------------------------------------------------------------    

! TEMPORARY
!      ELSE 
!        CALL RFLU_BLIN_ComputeXSectLine(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
!                                        iPatch,ifg,nXSect,xsd)        
!                              
!        IF ( nXSect > 0 ) THEN 
!          dist = xsd(1)
!          
!          IF ( dist < distOut ) THEN 
!            distOut = dist
!            iflOut  = ifl
!          END IF ! dist           
!        END IF ! nXSect         
!      END IF ! pPatch%bf2v        
! END TEMPORARY

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

    ELSE ! Defensive programming
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! iPatch
  END DO ! ifl

! ******************************************************************************
! Set output
! ******************************************************************************

  ipatchout = pc2f(1,iflout)
  ifgout    = pc2f(2,iflout)      

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_computelinecellxsectfast








! ******************************************************************************
!
! Purpose: Compute intersection of given line vector and faces of given cell
!   and distance between given location and intersection using safe algorithm. 
!
! 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
!   icg                 Global cell index
!
! Output:
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of intersection
!   distOut             Distance from location to intersection 
!   iPatchOut           Face location 
!   ifgOut              Face index 
!
! 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. Distance might be zero if vertices lie on path, so need to include this
!      case in IF statement on dist after calls to routine which computes the
!      distance.
!   3. This is called a safe algorithm because it can detect and correct 
!      inconsistent input data.    
!
! ******************************************************************************

SUBROUTINE rflu_computelinecellxsectsafe(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
                                         distout,ipatchout,ifgout)

  USE rflu_modbilinearpatch, ONLY: rflu_blin_computexsectline   
     
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: icg
  INTEGER, INTENT(OUT) :: ifgout,ipatchout
  REAL(RFREAL), INTENT(IN) :: ex,ey,ez
  REAL(RFREAL), INTENT(INOUT) :: xloc,yloc,zloc    
  REAL(RFREAL), INTENT(OUT) :: distout
  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: c1,c2,errorflag,icl,ict,ifl,iflout,iflout1,iflout2,ifg,ipatch, & 
             nfaces,nxsect
  INTEGER, DIMENSION(:,:), POINTER :: pc2f
  REAL(RFREAL) :: distout1,distout2,fnx,fny,fnz,toler,xcofg,ycofg,zcofg
  REAL(RFREAL), DIMENSION(2) :: xsd
  REAL(RFREAL), DIMENSION(6) :: denom,numer
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  global => pregion%global

  CALL registerfunction(global,'RFLU_ComputeLineCellXSectSafe',&
  'RFLU_ModGeometryTools.F90')

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

  pgrid => pregion%grid

  toler = -pregion%mixtInput%tolerICT ! Must be consistent with ICT tolerance

! ******************************************************************************
! 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)

! DEBUG
!  WRITE(0,*) '@@@100',pRegion%iRegionGlobal,'xsect with cell:',icg,ict,icl
! END DEBUG
                              
! ******************************************************************************
! Loop over faces of cell
! ******************************************************************************
  
  DO ifl = 1,nfaces
    ipatch = pc2f(1,ifl)
    ifg    = pc2f(2,ifl)

! DEBUG
!    WRITE(0,*) '@@@200',pRegion%iRegionGlobal,'Comp xsect with face:',iPatch,ifg
!        IF ( iPatch == 0 ) THEN
!          WRITE(0,*) '@@@210',pRegion%iRegionGlobal,pGrid%fn(1:3,ifg)
!        ELSE
!          pPatch => pRegion%patches(iPatch)
!          WRITE(0,*) '@@@211',pPatch%fn(1:3,ifg)         
!        END IF ! iPatch
! END DEBUG

! ==============================================================================
!   Interior face
! ==============================================================================

    IF ( ipatch == 0 ) THEN 
    
! ------------------------------------------------------------------------------
!     Triangular face
! ------------------------------------------------------------------------------    
    
! TEMPORARY
!      IF ( pGrid%f2v(4,ifg) == VERT_NONE ) THEN 
! END TEMPORARY
        c1 = pgrid%f2c(1,ifg)
        c2 = pgrid%f2c(2,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                    

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

        numer(ifl) = (xcofg-xloc)*fnx + (ycofg-yloc)*fny + (zcofg-zloc)*fnz
        denom(ifl) = ex*fnx + ey*fny + ez*fnz 

! ------------------------------------------------------------------------------
!     Quadrilateral face
! ------------------------------------------------------------------------------    
    
! TEMPORARY
!      ELSE 
!        CALL RFLU_BLIN_ComputeXSectLine(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
!                                        iPatch,ifg,nXSect,xsd)        
!
!! DEBUG
!        WRITE(0,*) '@@@300',pRegion%iRegionGlobal,'xsect dist:',nXSect,xsd
!! END DEBUG
!                              
!        IF ( nXSect > 0 ) THEN 
!          numer(ifl) = xsd(1)
!          denom(ifl) = 1.0_RFREAL
!        ELSE 
!          numer(ifl) = HUGE(1.0_RFREAL)
!          denom(ifl) = 1.0_RFREAL 
!        END IF ! nXSect  
!      END IF ! pGrid%f2v     
! END TEMPORARY

! ==============================================================================
!   Boundary face
! ==============================================================================

    ELSE IF ( ipatch > 0 ) THEN 
      ppatch => pregion%patches(ipatch)

! ------------------------------------------------------------------------------
!     Triangular face
! ------------------------------------------------------------------------------    
    
!      IF ( pPatch%bf2v(4,ifg) == VERT_NONE ) THEN 
        fnx = ppatch%fn(xcoord,ifg)
        fny = ppatch%fn(ycoord,ifg)
        fnz = ppatch%fn(zcoord,ifg)          

        xcofg = ppatch%fc(xcoord,ifg)
        ycofg = ppatch%fc(ycoord,ifg)
        zcofg = ppatch%fc(zcoord,ifg)
      
        numer(ifl) = (xcofg-xloc)*fnx + (ycofg-yloc)*fny + (zcofg-zloc)*fnz
        denom(ifl) = ex*fnx + ey*fny + ez*fnz      

! ------------------------------------------------------------------------------
!     Quadrilateral face
! ------------------------------------------------------------------------------    
    
! TEMPORARY
!      ELSE 
!        CALL RFLU_BLIN_ComputeXSectLine(pRegion,xLoc,yLoc,zLoc,ex,ey,ez,icg, & 
!                                        iPatch,ifg,nXSect,xsd)        
!                              
!        IF ( nXSect > 0 ) THEN 
!          numer(ifl) = xsd(1)
!          denom(ifl) = 1.0_RFREAL
!        ELSE 
!          numer(ifl) = HUGE(1.0_RFREAL)
!          denom(ifl) = 1.0_RFREAL 
!        END IF ! nXSect      
!      END IF ! pPatch%bf2v
! END TEMPORARY
    ELSE ! Defensive programming
      CALL errorstop(global,err_reached_default,__line__)
    END IF ! iPatch
  END DO ! ifl

! ******************************************************************************
! Set output
! ******************************************************************************

  errorflag = err_none

  distout1 = huge(1.0_rfreal)
  distout2 = huge(1.0_rfreal)

  iflout1 = crazy_value_int
  iflout2 = crazy_value_int

  faceloop: DO ifl = 1,nfaces
    IF ( numer(ifl) > toler ) THEN ! Inside cell icl
      IF ( denom(ifl) > 0.0 ) THEN
        distout = max(numer(ifl),0.0_rfreal)/denom(ifl)
        
        IF ( distout < distout1 ) THEN 
          distout2 = distout1        
          iflout2  = iflout1
          
          distout1 = distout
          iflout1  = ifl
        ELSE IF ( distout < distout2 ) THEN 
          distout2 = distout
          iflout2  = ifl        
        END IF ! dist
      END IF ! denom 
    ELSE ! Outside cell icl
      errorflag = 1
      distout = 0.0_rfreal
      iflout  = ifl
      
      EXIT faceloop
    END IF ! numer
  END DO faceloop

  IF ( errorflag == err_none ) THEN
    distout = distout1 ! Temporary 
    iflout  = iflout1     
  END IF ! errorFlag

! DEBUG
!  WRITE(0,*) '@@@400',pRegion%iRegionGlobal,iflOut,distOut
! END DEBUG
                              
  ipatchout = pc2f(1,iflout)
  ifgout    = pc2f(2,iflout)      

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_computelinecellxsectsafe






! ******************************************************************************
!
! Purpose: Determine whether point is inside bounding box.
!
! Description: 
!
! Input:
!   global              Pointer to global data
!   xLoc,yLoc,zLoc      x-, y-, and z-coordinates of location in question
!   xMin,yMin,zMin      Minimum x-, y-, and z-coordinates
!   xMax,yMax,zMax      Maximum x-, y-, and z-coordinates
!
! Output: 
!   RFLU_TestInBoundingBox = .TRUE.     if point inside bounding box
!   RFLU_TestInBoundingBox = .FALSE.    if point outside bounding box
!
! Notes: None.
!
! ******************************************************************************

LOGICAL FUNCTION rflu_testinboundbox(global,xLoc,yLoc,zLoc,xMin,xMax, & 
                                     ymin,ymax,zmin,zmax)
   
  IMPLICIT NONE

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

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

  REAL(RFREAL), INTENT(IN) :: xloc,xmax,xmin,yloc,ymax,ymin,zloc,zmax,zmin
  TYPE(t_global), POINTER :: global

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

  CALL registerfunction(global,'RFLU_TestInBoundBox',&
  'RFLU_ModGeometryTools.F90')

! ******************************************************************************
! Determine whether point is in bounding box
! ******************************************************************************

  rflu_testinboundbox = .true. 

  IF ( (xloc < xmin) .OR. (xloc > xmax) ) THEN 
    rflu_testinboundbox = .false. 
    CALL deregisterfunction(global)    
    RETURN 
  END IF ! xLoc

  IF ( (yloc < ymin) .OR. (yloc > ymax) ) THEN 
    rflu_testinboundbox = .false. 
    CALL deregisterfunction(global)      
    RETURN 
  END IF ! yLoc

  IF ( (zloc < zmin) .OR. (zloc > zmax) ) THEN 
    rflu_testinboundbox = .false. 
    CALL deregisterfunction(global)      
    RETURN 
  END IF ! zLoc

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

  CALL deregisterfunction(global)

END FUNCTION rflu_testinboundbox





! ******************************************************************************
!
! Purpose: Determine whether vector is aligned with Cartesian coordinate 
!   direction.
!
! Description: 
!
! Input:
!   global              Pointer to global data
!   dr                  Vector
!   dir                 Cartesian direction index
!
! Output: 
!   RFLU_TestVectorCartAxisAligned = .TRUE.     if dr aligned with Cartesian 
!                                               coordinate direction
!   RFLU_TestVectorCartAxisAligned = .FALSE.    if dr not aligned with Cartesian 
!                                               coordinate direction
!
! Notes: None.
!
! ******************************************************************************

LOGICAL FUNCTION rflu_testvectorcartaxisaligned(global,dr,dir)
   
  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: dir
  REAL(RFREAL), INTENT(INOUT) :: dr(xcoord:zcoord)
  TYPE(t_global), POINTER :: global

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

  REAL(RFREAL) :: drsum,drtol

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

  CALL registerfunction(global,'RFLU_TestVectorCartAxisAligned',&
  'RFLU_ModGeometryTools.F90')

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

  drtol = 1.0e-12_rfreal

! ******************************************************************************
! Test whether vector aligned 
! ******************************************************************************

  rflu_testvectorcartaxisaligned = .false. 

  dr(dir) = 0.0_rfreal
        
  drsum = abs(dr(xcoord)) + abs(dr(ycoord)) + abs(dr(zcoord))

  IF ( drsum <= drtol ) THEN 
    rflu_testvectorcartaxisaligned = .true.
  END IF ! drSum

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

  CALL deregisterfunction(global)

END FUNCTION rflu_testvectorcartaxisaligned





! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modgeometrytools


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModGeometryTools.F90,v $
! Revision 1.5  2008/12/06 08:44:22  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:33  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2007/03/06 18:06:42  haselbac
! Added function to test for alignment ofvector with Cartesian direction
!
! Revision 1.2  2006/04/07 15:19:19  haselbac
! Removed tabs
!
! Revision 1.1  2005/12/24 21:17:50  haselbac
! Initial revision
!
! ******************************************************************************