tem_determine_discreteVector Subroutine

public subroutine tem_determine_discreteVector(vector, compareVector, angle)

Compare the incoming discrete vector against a set of prevailing vectors and return the found closest prevailing integer vector

A set of real numbered unit vectors (compareVector) must be given. The integer numbered vector is compared against those and the closest integer vector overwrites the original vector.

Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: vector(:)

The given vector, will be set to the vector in compareVector that is best aligned to it.
real(kind=rk), intent(in) :: compareVector(:,:)

Set of unit vectors to select from against. size is (3, nVectors)
real(kind=rk), intent(out), optional :: angle

angle between vectors

Called by

proc~~tem_determine_discretevector~~CalledByGraph proc~tem_determine_discretevector tem_determine_discreteVector proc~tem_stencil_getlabelforcxdir tem_stencil_getLabelForcxDir proc~tem_stencil_getlabelforcxdir->proc~tem_determine_discretevector

Contents

Source Code

tem_determine_discreteVector

Source Code

  subroutine tem_determine_discreteVector( vector, compareVector, angle )
    ! -------------------------------------------------------------------- !
    !> The given vector, will be set to the vector in compareVector
    !! that is best aligned to it.
    integer, intent(inout) :: vector(:)
    !> Set of unit vectors to select from against.
    !! size is (3, nVectors)
    real( kind=rk ), intent(in) :: compareVector(:,:)
    !> angle between vectors
    real( kind=rk ), intent(out), optional :: angle
    ! -------------------------------------------------------------------- !
    integer :: iDir, bestIndex
    real(kind=rk) :: length, rv(3)
    real(kind=rk) :: rVectIn(3), angleLoc, dotProduct, maxplen, min_nz_comp
    ! -------------------------------------------------------------------- !

    angleLoc = 0.0_rk

    ! Only need to do anything for non-zero vectors.
    if ( maxval(abs(vector)) > 0 ) then
      rv = real(vector, kind=rk)

      length = sqrt(rv(1)**2 + rv(2)**2 + rv(3)**2)

      ! Normalize vector
      rVectIn = rv / length

      ! Keep track of the best fitting index.
      bestIndex = 0

      ! Keep track of the maximal projected length.
      ! As we are looking at two unit vectors, the smallest
      ! projected length should be -1 (opposing directions).
      ! Thus, starting with -2 here as a save comparison that is
      ! smaller than all possible outcomes from compareVectors.
      maxplen = -2.0_rk

      ! Go over all vectors in the compareVector
      do iDir = 1, size(compareVector, 2)
        ! Go through all prevailing directions and check the minimum
        dotProduct = rVectIn(1)*compareVector( 1, iDir ) &
          &        + rVectIn(2)*compareVector( 2, iDir ) &
          &        + rVectIn(3)*compareVector( 3, iDir )

        dotProduct = min( max(dotProduct, -1.0_rk), 1.0_rk )
        if ( dotProduct > maxplen ) then
          maxplen = dotProduct
          bestIndex = iDir
          if ( maxplen .feq. 1._rk ) EXIT
        end if
      end do

      ! Align the vector to the closest matching prevailing direction from
      ! the set of compareVector:

      ! Use the smallest nonzero component to scale the unit vector to
      ! components close to integer numbers.
      min_nz_comp = abs( minval( pack(compareVector(:, bestIndex),     &
        &                             abs(compareVector(:, bestIndex)) &
        &                             > epsilon(0.0_rk))               &
        &                      )                                       )
      ! Use the nearest integers to represent the projected vector.
      vector = nint(compareVector(:, bestIndex)/min_nz_comp)

      ! Angle between the original vector and the projected direction.
      angleLoc = acos(maxplen)
    end if ! not zero-vector

    if (present(angle)) angle = angleLoc

  end subroutine tem_determine_discreteVector