matrix_calculations.f90

Go to the documentation of this file.

module matrix_calculations

  use itm_types
  use mod_corners, only : nodeno
  use mod_helena_io

  implicit none

  contains

  subroutine set_node_number(n_node, n1, n2, n3, n4)

    implicit none

    integer(itm_i4), intent(in)  :: n_node
    integer(itm_i4), intent(out) :: n1, n2, n3, n4

    if (n_node .lt. lbound(nodeno,1)) then
      write(iu6,*) 'n_node out of bounds = ', n_node
      read(*,*)
      stop 'Error in set_node_number'
    else if (n_node .gt. ubound(nodeno,1)) then
      write(iu6,*) 'n_node out of bounds = ', n_node
      read(*,*)
      stop 'Error in set_node_number'
    end if

    n1 = nodeno(n_node, 1)
    n2 = nodeno(n_node, 2)
    n3 = nodeno(n_node, 3)
    n4 = nodeno(n_node, 4)

    return
  end subroutine set_node_number

  subroutine set_index(indices, inv)
!-----------------------------------------------------------------------
! Calculate the index in the matrix or its inverse. It is reordered from
! clockwise ordering in the rest of helena to help reduce the matrix size
! by a factor of 2.
! The poloidal index j changes to:
!          jn = 1                 (j = 1, np)
!          jn = 2 * (j - 1),      (2 .le. j .le. (np - 1) / 2 + 1)
!          jn = 2 * (np - j) + 1, ((np - 1) / 2 + 2 .le. j .le. np - 1)
!-----------------------------------------------------------------------
!
! the inverse of index in build_matrix will be restored
!-----------------------------------------------------------------------

    use mod_mesh, only : nr, np, ias

    implicit none

    integer(itm_i4), dimension(nr * np), intent(inout) :: indices
    logical, intent(in) :: inv
    integer(itm_i4) :: i, j, jn, ij1, ijn

    if (ias == 1) then
      do i = 1, nr
        j = 1
        jn = 1
        ij1 = (i - 1) * np + j
        ijn = (i - 1) * (np - 1) + jn
        if (inv) then
          indices(ijn) = ij1
        else
          indices(ij1) = ijn
        end if
        do j = 2, (np - 1) / 2 + 1
          jn = 2 * (j - 1)
          ij1 = (i - 1) * np + j
          ijn = (i - 1) * (np - 1) + jn
          if (inv) then
            indices(ijn) = ij1
          else
            indices(ij1) = ijn
          end if
        end do
        do j = (np - 1) / 2 + 2, np - 1
          jn = 2 * (np - j) + 1
          ij1 = (i - 1) * np + j
          ijn = (i - 1) * (np - 1) + jn
          if (inv) then
            indices(ijn) = ij1
          else           
            indices(ij1) = ijn
          end if
        end do

        if ( .not. inv )  then
          j = np
          jn = 2 * (np - j) + 1
          ij1 = (i - 1) * np + j
          ijn = (i - 1) * (np - 1) + jn
          indices(ij1) = ijn
        end if
      end do
      if (inv) return
    else
      do i = 1, nr
        do j = 1, np
          ij1 = (i - 1) * np + j
          indices(ij1) = ij1
        end do
      end do
    endif

    return
  end subroutine set_index

end module matrix_calculations