psmile_control_cell_2d_dble.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_Control_cell_2d_dble
!
! !INTERFACE:

      subroutine psmile_control_cell_2d_dble (grid_id,               &
                      ic, nc, icdim1, list, j, k,                    &
                      coords1,  coords2,  shape,                     &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, fnd, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Control_cell_2d_dble
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id
!
!     Grid Id (for info degraded cells)
!
      Integer, Intent (In)            :: nc
!
!     Number of cells (and points) to be controlled
!
      Integer, Intent (In)            :: icdim1
!
!     First dimension of "ic"
!
      Integer, Intent (In)            :: ic (icdim1, ndim_2d)
!
!     Indices of the (corner) cells to be controlled
!     ic (1:nc, 1) = 1st index (i-index) of the corner cell
!     ic (1:nc, 2) = 2nd index (j-index) of the corner cell
!
      Integer, Intent (In)            :: list (nc)
!
!     I-Indices of the points in to be controlled
!
      Integer, Intent (In)            :: j
!
!     J-Index of the points in to be controlled
!
      Integer, Intent (In)            :: k
!
!     K-Index of the points in to be controlled
!
      Integer, Intent (In)            :: shape (2, ndim_3d)

!     Dimension of coordinate arrays "coord1" and "coord2"
!     which contain the coordinates to be searched.

      Double Precision, Intent (In)   :: coords1 (shape(1,1):shape(2,1), 
                                                  shape(1,2):shape(2,2), 
                                                  shape(1,3):shape(2,3))
      Double Precision, Intent (In)   :: coords2 (shape(1,1):shape(2,1), 
                                                  shape(1,2):shape(2,2), 
                                                  shape(1,3):shape(2,3))

!     Coordinates to be searched

      Integer, Intent (In)            :: corner_shape (2, ndim_2d)

!     Dimension of corner coordinate arrays

      Integer, Intent (In)            :: nbr_corners

!     Number of corners
!     Note: This routine assumes nbr_corners == 4.

      Double Precision,   Intent (In) ::                 
         corners1 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

      Double Precision,   Intent (In) ::                 
         corners2 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

!     Corner coordinates

      Double Precision, Intent (In)   :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0
!
! !OUTPUT PARAMETERS:
!
      Logical, Intent (Out)           :: fnd (nc)
!
!     Was the i-th point of "list" found within cell "ic"
!
      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILe_Control_cell_2d_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NTET = Number of triangles used to determine whether
!         point is within the grid cell
!  EPS  = Tolerance for determinate of deformed grid cells
!
      Integer, Parameter              :: ntet = 3
!
      Double Precision, Parameter     :: eps = 1.0d-9
      Double Precision, Parameter     :: rmxeps = 1.0d30
      Double Precision, Parameter     :: dzero = 0.0
!
! !LOCAL VARIABLES
!
!     ... for tolerances
!
      Double Precision                :: tol1
!
!     ... for locations controlled
!
      Integer                         :: i, jj, m, n
!
!     ... Logical reference point and edges (constant)
! 
      Integer, Save                   :: ijkedg (2, ntet)
!
!     ... Reference point and edges and distances
! 
      Double Precision                :: xyzpnt (nc, ntet, ndim_2d), 
                                         xyzref (nc, ntet, ndim_2d)
      Double Precision                :: xyzedg (nc, ntet, ndim_2d, 2)
      Double Precision                :: d1, d2
      Double Precision                :: det (nc, ntet)
      Double Precision                :: dnorm
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Control_cell_2d_dble" controls whether the point
! "(list(i),j,k)" is contained in the corner cell "ic(i,*)".
! The routine assumes that the number of corners for the 2-dimensional
! cell is 4.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 12.03.07    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_control_cell_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_control_cell_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative indices of triangles
!
!     The cell is split into 3 triangles
!
!                   4-------3
!                   !       !
!                   !       !
!            Corner 1 ----- 2
!
      data ((ijkedg (jj, n), jj = 1, 2), n = 1, ntet) &
           / 2, 4,  3, 1,  4, 2/
!   Corner    1      2      3
!   Edges of corner 1 are connections to Corner 2 and Corner 4
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), nc

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (tol < 0.0) then
         call psmile_assert (__FILE__, __LINE__, &
                             "tol must be >= 0.0")
      endif
!
!cdir vector
         do i = 1, nc
         if (ic(i,1) < corner_shape(1,1) .or. ic(i,1) > corner_shape(2,1) .or. &
             ic(i,2) < corner_shape(1,2) .or. ic(i,2) > corner_shape(2,2) ) then
            exit
         endif
         end do
!
      if (i <= nc) then
         print *, 'i, ic', i, ic(i, :)
         print *, 'corner_shape', corner_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index')
      endif
#endif
!
      ierror = 0
!
!===> Compute reference points, corners and edges of transformed
!     triangles.
!
        do n = 1, ntet
!cdir vector
           do i = 1, nc
           xyzref (i, n, 1) = corners1 (ic(i,1), ic(i,2), n)
           xyzref (i, n, 2) = corners2 (ic(i,1), ic(i,2), n)
           end do
!
        xyzpnt (:, n, 1) = coords1 (list(:), j, k) - xyzref (:, n, 1)
        xyzpnt (:, n, 2) = coords2 (list(:), j, k) - xyzref (:, n, 2)
!
           do m = 1, 2
!cdir vector
              do i = 1, nc
              xyzedg (i, n, 1, m) = corners1 (ic(i,1), ic(i,2), ijkedg(m, n)) &
                                  - xyzref (i, n, 1)

              xyzedg (i, n, 2, m) = corners2 (ic(i,1), ic(i,2), ijkedg(m, n)) &
                                  - xyzref (i, n, 2)
              end do
           end do
        end do ! n
!
!===> Solve bilinear system; compute spat product
!
!   (xyzedg(*,1), xyzedg(*,2)) (x) = xyzpnt
!
      det(:, :) = xyzedg(:, :, 1,1)*xyzedg(:, :, 2,2) &
                - xyzedg(:, :, 1,2)*xyzedg(:, :, 2,1)

#ifdef DEBUG_TRACE
      if (j == ictl_ind(2)) then
!cdir vector
            do i = 1, nc
            if ( list(i) == ictl_ind(1) ) exit
            end do

         if (i <= nc) then
            print *, 'ictl: i,j', list(i), j, coords1 (list(i), j, k), coords2 (list(i), j, k)
            do n = 1, 4
               print *, corners1 (ic(i,1), ic(i,2), n), corners2 (ic(i,1), ic(i,2), n)
            enddo
            print *, '#################################################'
            print *, ' i, nc ', i, nc
            print *, ' det', det (i, :)
         end if
      end if
#endif
!
#define ALLOW_DEGRADED_CELL
#ifdef ALLOW_DEGRADED_CELL
!
! ??? ist das koscher ?? wird d1(n), d2(n), d3(n) == 0
!
         do n = 1, ntet
!cdir vector
            do i = 1, nc
            dnorm = abs(xyzedg(i, n,1,1)) + abs(xyzedg(i, n,2,1)) + &
                    abs(xyzedg(i, n,1,2)) + abs(xyzedg(i, n,2,2))

            if (abs(det(i, n)) < eps*dnorm*dnorm .or. dnorm == dzero) then
               det (i, n) = sign (rmxeps, det(i, n))
            else
               det (i, n) = 1.0d0 / det (i, n)
            endif
            end do ! i
         end do ! n

#else 

loop_n:  do n = 1, ntet
!cdir vector
            do i = 1, nc
            dnorm = abs(xyzedg(i, n,1,1)) + abs(xyzedg(i, n,2,1)) + &
                    abs(xyzedg(i, n,1,2)) + abs(xyzedg(i, n,2,2))

            if (abs(det(i, n)) < eps*dnorm*dnorm .or. &
                dnorm(i,n) == dzero) then
               det (i, n) = dnorm * 1e-25
! eigentlich exit loop_n und print_degraded cell
! ist bei vektorisierung verschwunden
! sollte man in part_list suchen
            endif
            end do ! i
         end do ! n
!
#ifdef PRINT_DEGRADED_CELL
      if (n < ntet) then
!
!===> cell is degraded
!
         print 9960, trim(ch_id), grid_id, ic (i, :), &
               corners1(ic(i,1), ic(i,2), 1),  &
               corners2(ic(i,1), ic(i,2), 1),  &
               corners1(ic(i,1), ic(i,2), 2),  &
               corners2(ic(i,1), ic(i,2), 2),  &
               corners1(ic(i,1), ic(i,2), 3),  &
               corners2(ic(i,1), ic(i,2), 3),  &
               corners1(ic(i,1), ic(i,2), 4),  &
               corners2(ic(i,1), ic(i,2), 4)
!
            do n = 1, ntet
            if (abs(det(i, n)) < eps) then
               print 9950, n, det (i, n), dnorm, &
                 xyzedg (i, n,1, 1), xyzedg (i, n,1, 2), &
                 xyzedg (i, n,2, 1), xyzedg (i, n,2, 2)
            endif
            end do
!
         nicht gefunden fnd (i) = .false. versucht, durch
         ganz kleine Determinante zu loesen
      endif
#endif /* PRINT_DEGRADED_CELL */
!
      det (:, :) = 1.0d0 / det (:, :)
!
#endif /* ALLOW_DEGRADED_CELL */
!
!===> Return mask for points "list(i)"
!     which were found in the specific cell (ic(i,*), j,k)
!
      tol1 = tol + 1.0
      fnd (:) = .false.
!
         do n = 1, ntet
!cdir vector
            do i = 1, nc

            d1 = (xyzpnt(i, n,1)*xyzedg(i, n,2,2)- &
                  xyzpnt(i, n,2)*xyzedg(i, n,1,2)) * det (i, n)
            d2 = (xyzpnt(i, n,2)*xyzedg(i, n,1,1)- &
                  xyzpnt(i, n,1)*xyzedg(i, n,2,1)) * det (i, n)
!
            fnd (i) = fnd (i) .or. &
               (min(d1, d2) >= -tol .and. &
                max(d1,dzero) + max(d2,dzero) <= tol1)
            end do ! nc
         end do ! ntet
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), count(fnd(:))

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

9990 format (1x, a, ': psmile_control_cell_2d_dble: nc', i5)
9980 format (1x, a, ': psmile_control_cell_2d_dble: eof, suc', i5)

#ifdef PRINT_DEGRADED_CELL
9960 format (1x, a, ': psmile_control_cell_2d_dble: grid id =', i5, &
            /1x,  ' Grid cell (', i4, ',', i4,                    &
                 ') is degraded! Cell coordinates are:', 1p,      &
            4 (/1x, '(', 1 (e17.9, ','), e17.9, ')'))
9950 format (1x, i1, '-th spawn; det =', 1p, e17.9, ' (', e17.9, &
             ') of', 2 (/1x, 2e17.9, ','))
#endif

      end subroutine psmile_control_cell_2d_dble